🧪 Historia powstania zbioru danych QM9¶

Zbiór danych QM9 jest jednym z najważniejszych i najczęściej wykorzystywanych benchmarków w chemoinformatyce, uczeniu maszynowym dla chemii oraz projektowaniu cząsteczek. Jego stworzenie było odpowiedzią na potrzebę posiadania dużego, jednorodnego i wysokiej jakości zbioru danych kwantowo-chemicznych, który umożliwiłby rozwój algorytmów ML w obszarze właściwości molekuł.

---

📜 Geneza projektu¶

Zbiór QM9 został opracowany w latach 2014–2017 przez zespół naukowców z University of Basel, Cornell University i University of Vienna: Ramakrishnan, Dral, Rupp i von Lilienfeld. Został po raz pierwszy opisany w publikacji:

“Quantum chemistry structures and properties of 134 kilo molecules” (Scientific Data, 2014)

Celem autorów było stworzenie systematycznego i kompletnego zbioru obejmującego tysiące małych cząsteczek organicznych oraz ich właściwości obliczonych metodami chemii kwantowej.

---

🧩 Jak powstał zbiór?¶

Proces tworzenia QM9 obejmował kilka etapów:

1. Wybór cząsteczek (enumeracja strukturalna)¶

Autorzy wybrali wszystkie stabilne cząsteczki organiczne o wzorze sumarycznym zawierającym tylko:

• C,
• H,
• O,
• N,
• F, oraz maksymalnie 9 atomów ciężkich. Dlatego nazwa datasetu to QM9.

📌 W sumie zebrano 133 885 struktur.

---

2. Optymalizacja geometrii¶

Dla każdej cząsteczki wykonano obliczenia metodą DFT (Density Functional Theory) z funkcjonałem:

• B3LYP
• bazą 6-31G(2df,p)

Była to ogromna kampania obliczeniowa — w sumie miliony iteracji obliczeń kwantowych.

---

3. Wyznaczenie właściwości fizykochemicznych¶

Dla każdej cząsteczki obliczono ≥ 17 właściwości, m.in.:

• energie atomizacji,
• moment dipolowy,
• HOMO, LUMO, gap,
• częstotliwości drgań,
• polaryzowalność,
• masę,
• geometrię 3D,
• temperaturę wrzenia (przybliżoną metodami kwantowymi).

📊 Dzięki temu zbiór jest idealny do regresji, predykcji właściwości i zadań ML w chemii.

---

4. Standaryzacja i walidacja¶

Każdą strukturę zweryfikowano:

• poprawność SMILES,
• poprawność minimalnej energii,
• brak fragmentacji,
• poprawna topologia.

Wczesne wersje datasetu zawierały ~300 błędnych geometrii — zostały później usunięte.

---

🚀 Znaczenie QM9 dla współczesnego ML w chemii¶

Zbiór QM9 stał się złotym standardem w testowaniu modeli:

• GNN (Graph Neural Networks),
• Message Passing Neural Networks,
• SchNet, PhysNet, DimeNet,
• Transformerów molekularnych,
• modeli generatywnych (VAE, diffusion models).

Ponieważ dane są jednorodne, czyste i precyzyjne, QM9 pozwala porównywać algorytmy w sposób uczciwy i powtarzalny.

---

⭐ Dlaczego QM9 jest tak ważny?¶

• zawiera największy do tamtej pory publiczny zbiór obliczeń DFT,
• jest kompletny w zakresie małych cząsteczek organicznych,
• ma wysoką jakość,
• idealnie nadaje się do projektów ML w chemii ❤️‍🔥.

Dziś stanowi fundament ogromnej części badań nad AI for Molecules.

---

📂 Opis zbioru danych QM9¶

W tej sekcji przedstawiono szczegółowy opis struktur danych dostępnych w roboczej wersji zbioru QM9, zawierającej 1000 losowo wybranych cząsteczek. Dane pochodzą z pełnego zbioru 133 885 struktur obliczonych metodami DFT.

---

🧬 Struktura tabeli danych¶

Po wczytaniu pliku dataset zawiera następujące typowe kolumny (nazwy mogą różnić się zależnie od wersji):

Kolumna	Opis	Typ
smiles	SMILES cząsteczki (struktura 2D zakodowana tekstowo)	string
formula	Wzór sumaryczny cząsteczki	string
rotational_constants_A/B/C	Stałe rotacyjne (GHz) wyliczone z geometrii 3D	float
dipole_moment	Moment dipolowy (Debye)	float
isotropic_polarizability	Izotropowa polaryzowalność (Bohr³)	float
homo	Energia HOMO (eV)	float
lumo	Energia LUMO (eV)	float
gap	Przerwa energetyczna HOMO–LUMO (eV)	float
electronic_spatial_extent	Objętość elektroniczna ⟨r²⟩ (Bohr²)	float
zero_point_energy	Energia zerowa (ZPE, Hartree)	float
internal_energy_0K	Energia wewnętrzna w 0 K (Hartree)	float
internal_energy_298K	Energia w 298.15 K (Hartree)	float
enthalpy_298K	Entalpia w 298.15 K (Hartree)	float
free_energy_298K	Energia swobodna (Gibbs) w 298.15 K (Hartree)	float
heat_capacity	Ciepło właściwe Cv (cal/mol·K)	float
mulliken_charges (opcjonalne)	Rozmieszczenie ładunku Mullikena	float/list
coordinates_X/Y/Z (opcjonalne)	Współrzędne 3D atomów (Å)	list
atom_types	Lista atomów (C/H/O/N/F)	list
num_atoms	Liczba atomów	int
num_heavy_atoms	Liczba atomów ciężkich (C/O/N/F)	int
molecule_id	Unikatowy identyfikator cząsteczki	int

⚠️ Uwaga: Dokładne kolumny w Twoim Excelu zależą od wersji konwertowanego pliku QM9, ale powyższa lista obejmuje cały standardowy zakres używany w chemoinformatyce.

---

🧠 Kategorie cech w zbiorze¶

Aby lepiej zrozumieć strukturę danych, można podzielić kolumny QM9 na 4 główne grupy:

1️⃣ Cechy strukturalne¶

• SMILES
• Wzór sumaryczny
• Liczba atomów (ogółem i ciężkich)
• Współrzędne 3D
• Typy atomów

👉 używane głównie do modelowania GNN i obliczeń geometrycznych.

---

2️⃣ Cechy energetyczne (DFT)¶

• ZPE
• Energie wewnętrzne
• Enthalpia
• Energia swobodna
• HOMO / LUMO / GAP

🔥 kluczowe dla modeli predykcyjnych właściwości kwantowych.

---

3️⃣ Cechy spektralne i rotacyjne¶

• Stałe rotacyjne A/B/C
• Ciepło właściwe
• Przestrzenne rozciągnięcie ładunku (⟨r²⟩)

---

4️⃣ Cechy elektrostatyczne¶

• Moment dipolowy
• Polaryzowalność
• Potencjały elektronowe i ładunki (jeśli dostępne)

---

🎯 Dlaczego te dane są wartościowe?¶

• Pochodzą z precyzyjnych obliczeń kwantowych, nie eksperymentów.
• Charakteryzują się wysoką jednorodnością metodologiczną.
• Idealnie nadają się do trenowania modeli ML predykcji właściwości molekularnych.
• Zawierają różnorodność strukturalną przy ograniczonej wielkości cząsteczek — kluczowe do nauki modeli GNN.

🎯 Cel projektu: Eksploracyjna analiza danych (EDA) zbioru QM9¶

Celem tego projektu jest przeprowadzenie kompleksowej eksploracyjnej analizy danych (EDA) na zbiorze QM9, zawierającym struktury i właściwości kwantowo-chemiczne małych cząsteczek organicznych.

---

🔹 Główne założenia projektu¶

1. Zrozumienie struktury danych

   • Identyfikacja kolumn, typów danych, braków danych i wyjątków.
   • Analiza podziału danych na cechy strukturalne, energetyczne, elektrostatyczne i rotacyjne.

2. Analiza statystyczna właściwości molekuł

   • Rozkłady podstawowych parametrów: energie, moment dipolowy, HOMO/LUMO, GAP.
   • Wykrywanie wartości odstających i potencjalnych błędów w danych.

3. Wizualizacja danych chemicznych

   • Graficzne przedstawienie relacji między właściwościami.
   • Wizualizacja cząsteczek na podstawie SMILES i geometrii 3D (jeśli środowisko pozwoli).

4. Tworzenie podzbiorów roboczych

   • Wydzielenie mniejszych próbek danych (np. 1000 wierszy) w celu szybkiej pracy i testowania modeli.

5. Przygotowanie danych do modeli ML

   • Identyfikacja cech predykcyjnych i docelowych.
   • Wstępne przekształcenia danych i normalizacja dla przyszłych zadań predykcyjnych.

---

🌟 Wartość projektu¶

• Projekt pozwala pokazać umiejętności w EDA i chemoinformatyce, w tym analizę dużych zbiorów danych chemicznych.
• Prezentuje profesjonalny workflow od wczytania danych, przez ich eksplorację, aż po wizualizację i przygotowanie do ML.
• Stanowi doskonały element portfolio dla kandydatów na stanowiska Data Scientist / Chemoinformatyk / Machine Learning Engineer w chemii.

In [3]:

import pandas as pd

# Ścieżka do dużego pliku
data = r"C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_dataset.xlsx"

# Wczytujemy z opcją ograniczenia do wybranych kolumn lub bez – jak wolisz
df = pd.read_excel(data)

# Losowa próbka 1000 wierszy z zachowaniem miksu danych
sample_df = df.sample(n=1000, random_state=42)

# Zapis nowego pliku roboczego
output_path = r"C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_sample_1000.xlsx"
sample_df.to_excel(output_path, index=False)

print("✔️ Utworzono plik roboczy z 1000 wierszy:")
print(output_path)

✔️ Utworzono plik roboczy z 1000 wierszy:
C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_sample_1000.xlsx

duzo recordów tworzenie roboczego df z 1000 wierszami¶

In [2]:

import pandas as pd

data = r"C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_sample_1000.xlsx"

df = pd.read_excel(data)
print(df.head())

                                                   X        y1        y2  \
0  <rdkit.Chem.rdchem.Mol object at 0x000001BF0ED... -1.291747  1.802935   
1  <rdkit.Chem.rdchem.Mol object at 0x000001BEFC9... -0.644026 -0.218808   
2  <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0... -0.459969  0.297628   
3  <rdkit.Chem.rdchem.Mol object at 0x000001BF2AB... -0.775036  1.089651   
4  <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0...  0.346542 -0.820930   

         y3        y4        y5        y6        y7        y8        y9  ...  \
0 -0.427511  0.711080  0.924869  1.096639  1.652419  1.864530  1.112639  ...   
1  2.287883  0.424803 -0.721806 -0.398991 -0.414869  0.044067 -0.463779  ...   
2  0.167480  0.130518  0.046509 -0.251763  0.650823  0.217649 -0.113869  ...   
3 -0.104059  1.107464  1.156964  0.315530  1.703451  1.364106  0.237038  ...   
4 -0.902704  0.618991  1.072929  0.839014 -0.128690  0.650888 -1.836564  ...   

   w5  w6  w7  w8  w9  w10  w11  w12  \
0   1   1   1   1   1    1    1    1   
1   1   1   1   1   1    1    1    1   
2   1   1   1   1   1    1    1    1   
3   1   1   1   1   1    1    1    1   
4   1   1   1   1   1    1    1    1   

                                                 ids  split  
0  [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C...  train  
1  [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@...  train  
2  [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]...  train  
3    [H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H]   test  
4  [H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2(...  train  

[5 rows x 27 columns]

zmiana nazw kolumn¶

dodaaćlink do data setu + usuniecie kolumn z wagami

In [3]:

# ============================================================
# KOMÓRKA 2: ZMIANA NAZW KOLUMN NA OPISOWE
# ============================================================

# Mapowanie oryginalnych nazw na opisowe
column_mapping = {
    # Identyfikatory i features
    'ids': 'SMILES',
    'X': 'SMILES_raw',

    # Właściwości molekularne (y1-y12)
    'y1': 'dipole_moment',          # μ (Debye) - moment dipolowy
    'y2': 'polarizability',         # α (Bohr³) - polaryzowalność
    'y3': 'HOMO',                   # εHOMO (eV) - energia HOMO
    'y4': 'LUMO',                   # εLUMO (eV) - energia LUMO
    'y5': 'gap',                    # Δε (eV) - HOMO-LUMO gap
    'y6': 'electronic_spatial_extent', # ⟨R²⟩ (Bohr²) - elektroniczny zasięg przestrzenny
    'y7': 'zpve',                   # ZPVE (eV) - zero point vibrational energy
    'y8': 'U0',                     # U₀ (eV) - energia wewnętrzna w 0K
    'y9': 'U',                      # U (eV) - energia wewnętrzna w 298.15K
    'y10': 'H',                     # H (eV) - entalpia w 298.15K
    'y11': 'G',                     # G (eV) - energia swobodna Gibbsa w 298.15K
    'y12': 'Cv',                    # Cv (cal/mol·K) - pojemność cieplna w 298.15K

    # Wagi (w1-w12)
    'w1': 'weight_dipole',
    'w2': 'weight_polarizability',
    'w3': 'weight_HOMO',
    'w4': 'weight_LUMO',
    'w5': 'weight_gap',
    'w6': 'weight_R2',
    'w7': 'weight_zpve',
    'w8': 'weight_U0',
    'w9': 'weight_U',
    'w10': 'weight_H',
    'w11': 'weight_G',
    'w12': 'weight_Cv',

    # Split (jeśli istnieje)
    'split': 'dataset_split'
}

# Zmiana nazw kolumn
df = df.rename(columns=column_mapping)

# ============================================================
# USUNIĘCIE KOLUMN Z WAGAMI (wszystkie = 1)
# ============================================================

# Lista kolumn z wagami do usunięcia
weight_columns = [col for col in df.columns if col.startswith('weight_')]

print("="*70)
print("USUWANIE KOLUMN Z WAGAMI")
print("="*70)
print(f"\n🗑️  Usuwam {len(weight_columns)} kolumn z wagami:")
for col in weight_columns:
    print(f"  • {col}")

# Usunięcie kolumn
df = df.drop(columns=weight_columns)

print(f"\n✅ Kolumny usunięte! Pozostało {len(df.columns)} kolumn")

print("\n" + "="*70)
print("ZMIENIONO NAZWY KOLUMN")
print("="*70)

print("\n📋 Nowe nazwy kolumn:")
for i, col in enumerate(df.columns, 1):
    print(f"  {i:2d}. {col}")

print("\n🔍 Pierwsze 3 wiersze po zmianie:")
print(df.head(3))

print("\n📊 Kształt DataFrame:")
print(f"  Wiersze: {len(df):,}")
print(f"  Kolumny: {len(df.columns)}")

print("\n✅ Nazwy kolumn zostały zmienione!")
print("="*70)

# Opcjonalnie: Wyświetl info o właściwościach
print("\n📖 OPIS WŁAŚCIWOŚCI KWANTOWO-CHEMICZNYCH:")
print("-" * 70)
properties_description = {
    'dipole_moment': 'Moment dipolowy [Debye] - miara polarności cząsteczki',
    'polarizability': 'Polaryzowalność [Bohr³] - odpowiedź na pole elektryczne',
    'HOMO': 'Energia HOMO [eV] - najwyższy zajęty orbital molekularny',
    'LUMO': 'Energia LUMO [eV] - najniższy niezajęty orbital molekularny',
    'gap': 'HOMO-LUMO gap [eV] - różnica energii, kluczowa dla reaktywności',
    'electronic_spatial_extent': 'Elektroniczny zasięg przestrzenny [Bohr²]',
    'zpve': 'Energia drgań punktu zerowego [eV]',
    'U0': 'Energia wewnętrzna w 0K [eV]',
    'U': 'Energia wewnętrzna w 298.15K [eV]',
    'H': 'Entalpia w 298.15K [eV]',
    'G': 'Energia swobodna Gibbsa w 298.15K [eV]',
    'Cv': 'Pojemność cieplna w 298.15K [cal/mol·K]'
}

for prop, desc in properties_description.items():
    print(f"  • {prop:30s} - {desc}")

print("="*70)

======================================================================
USUWANIE KOLUMN Z WAGAMI
======================================================================

🗑️  Usuwam 12 kolumn z wagami:
  • weight_dipole
  • weight_polarizability
  • weight_HOMO
  • weight_LUMO
  • weight_gap
  • weight_R2
  • weight_zpve
  • weight_U0
  • weight_U
  • weight_H
  • weight_G
  • weight_Cv

✅ Kolumny usunięte! Pozostało 15 kolumn

======================================================================
ZMIENIONO NAZWY KOLUMN
======================================================================

📋 Nowe nazwy kolumn:
   1. SMILES_raw
   2. dipole_moment
   3. polarizability
   4. HOMO
   5. LUMO
   6. gap
   7. electronic_spatial_extent
   8. zpve
   9. U0
  10. U
  11. H
  12. G
  13. Cv
  14. SMILES
  15. dataset_split

🔍 Pierwsze 3 wiersze po zmianie:
                                          SMILES_raw  dipole_moment  \
0  <rdkit.Chem.rdchem.Mol object at 0x000001BF0ED...      -1.291747   
1  <rdkit.Chem.rdchem.Mol object at 0x000001BEFC9...      -0.644026   
2  <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0...      -0.459969   

   polarizability      HOMO      LUMO       gap  electronic_spatial_extent  \
0        1.802935 -0.427511  0.711080  0.924869                   1.096639   
1       -0.218808  2.287883  0.424803 -0.721806                  -0.398991   
2        0.297628  0.167480  0.130518  0.046509                  -0.251763   

       zpve        U0         U         H         G        Cv  \
0  1.652419  1.864530  1.112639  1.112698  1.112698  1.112566   
1 -0.414869  0.044067 -0.463779 -0.463790 -0.463790 -0.463751   
2  0.650823  0.217649 -0.113869 -0.113867 -0.113867 -0.113876   

                                              SMILES dataset_split  
0  [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C...         train  
1  [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@...         train  
2  [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]...         train  

📊 Kształt DataFrame:
  Wiersze: 1,000
  Kolumny: 15

✅ Nazwy kolumn zostały zmienione!
======================================================================

📖 OPIS WŁAŚCIWOŚCI KWANTOWO-CHEMICZNYCH:
----------------------------------------------------------------------
  • dipole_moment                  - Moment dipolowy [Debye] - miara polarności cząsteczki
  • polarizability                 - Polaryzowalność [Bohr³] - odpowiedź na pole elektryczne
  • HOMO                           - Energia HOMO [eV] - najwyższy zajęty orbital molekularny
  • LUMO                           - Energia LUMO [eV] - najniższy niezajęty orbital molekularny
  • gap                            - HOMO-LUMO gap [eV] - różnica energii, kluczowa dla reaktywności
  • electronic_spatial_extent      - Elektroniczny zasięg przestrzenny [Bohr²]
  • zpve                           - Energia drgań punktu zerowego [eV]
  • U0                             - Energia wewnętrzna w 0K [eV]
  • U                              - Energia wewnętrzna w 298.15K [eV]
  • H                              - Entalpia w 298.15K [eV]
  • G                              - Energia swobodna Gibbsa w 298.15K [eV]
  • Cv                             - Pojemność cieplna w 298.15K [cal/mol·K]
======================================================================

In [4]:

df.head(10)

Out[4]:

	SMILES_raw	dipole_moment	polarizability	HOMO	LUMO	gap	electronic_spatial_extent	zpve	U0	U	H	G	Cv	SMILES	dataset_split
0	<rdkit.Chem.rdchem.Mol object at 0x000001BF0ED...	-1.291747	1.802935	-0.427511	0.711080	0.924869	1.096639	1.652419	1.864530	1.112639	1.112698	1.112698	1.112566	[H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C...	train
1	<rdkit.Chem.rdchem.Mol object at 0x000001BEFC9...	-0.644026	-0.218808	2.287883	0.424803	-0.721806	-0.398991	-0.414869	0.044067	-0.463779	-0.463790	-0.463790	-0.463751	[H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@...	train
2	<rdkit.Chem.rdchem.Mol object at 0x000001BF0E0...	-0.459969	0.297628	0.167480	0.130518	0.046509	-0.251763	0.650823	0.217649	-0.113869	-0.113867	-0.113867	-0.113876	[H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]...	train
3	<rdkit.Chem.rdchem.Mol object at 0x000001BF2AB...	-0.775036	1.089651	-0.104059	1.107464	1.156964	0.315530	1.703451	1.364106	0.237038	0.237068	0.237068	0.237016	[H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H]	test
4	<rdkit.Chem.rdchem.Mol object at 0x000001BF0E0...	0.346542	-0.820930	-0.902704	0.618991	1.072929	0.839014	-0.128690	0.650888	-1.836564	-1.836554	-1.836554	-1.836579	[H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2(...	train
5	<rdkit.Chem.rdchem.Mol object at 0x000001BF2A8...	0.409150	-2.145600	-0.758948	-0.494088	-0.113557	-1.271677	-1.391723	-1.357742	-0.503130	-0.503165	-0.503165	-0.503080	[H][N-]C1OC([H])([H])[C@]12[N@@H+]1C([H])([H])...	valid
6	<rdkit.Chem.rdchem.Mol object at 0x000001BF2AB...	-0.469601	0.294155	0.754484	0.791158	0.412659	-0.396482	0.637064	0.321416	-0.112035	-0.112037	-0.112037	-0.112018	[H]N([H])C1[NH2+][C@]2([H])C([H])([H])[C@]2(C(...	test
7	<rdkit.Chem.rdchem.Mol object at 0x000001BEFC5...	0.327348	-0.128489	-0.611199	-1.142716	-0.835853	-0.335116	-1.578459	-1.002448	-0.027435	-0.027466	-0.027466	-0.027395	[H]C1([H])[N@@H+]2[C@@]3([H])[C@]2([H])[C@@]2(...	train
8	<rdkit.Chem.rdchem.Mol object at 0x000001BF2A2...	1.261589	-2.541611	-1.797187	-0.974553	-0.073540	-1.369651	-1.949059	-2.116866	0.891255	0.891221	0.891221	0.891306	[H]C([H])([H])O[C@@]1(C([H])([H])[H])C([H])([H...	train
9	<rdkit.Chem.rdchem.Mol object at 0x000001BF0E1...	1.653432	-0.118068	0.462979	-1.759313	-1.990326	-0.729647	-1.155817	-1.150208	-0.435730	-0.435769	-0.435769	-0.435686	[H]C(=O)c1c([H])c(N([H])[H])nn1[H]	train

zmana 1 kolumny na licaby¶

In [5]:

# ============================================================
# KOMÓRKA 3: USUNIĘCIE ZBĘDNEJ KOLUMNY SMILES_raw
# ============================================================

print("="*70)
print("USUWANIE ZBĘDNEJ KOLUMNY")
print("="*70)

# Sprawdzenie czy kolumna istnieje
if 'SMILES_raw' in df.columns:
    print("\n🗑️  Usuwam kolumnę: SMILES_raw (duplikat SMILES)")
    df = df.drop(columns=['SMILES_raw'])
    print("✅ Kolumna usunięta!")
else:
    print("\n⚠️  Kolumna SMILES_raw nie istnieje")

print(f"\n📊 Aktualne kolumny ({len(df.columns)}):")
for i, col in enumerate(df.columns, 1):
    print(f"  {i:2d}. {col}")

print(f"\n📈 Kształt DataFrame: {df.shape}")
print("="*70)

======================================================================
USUWANIE ZBĘDNEJ KOLUMNY
======================================================================

🗑️  Usuwam kolumnę: SMILES_raw (duplikat SMILES)
✅ Kolumna usunięta!

📊 Aktualne kolumny (14):
   1. dipole_moment
   2. polarizability
   3. HOMO
   4. LUMO
   5. gap
   6. electronic_spatial_extent
   7. zpve
   8. U0
   9. U
  10. H
  11. G
  12. Cv
  13. SMILES
  14. dataset_split

📈 Kształt DataFrame: (1000, 14)
======================================================================

In [6]:

df.head(10)

Out[6]:

	dipole_moment	polarizability	HOMO	LUMO	gap	electronic_spatial_extent	zpve	U0	U	H	G	Cv	SMILES	dataset_split
0	-1.291747	1.802935	-0.427511	0.711080	0.924869	1.096639	1.652419	1.864530	1.112639	1.112698	1.112698	1.112566	[H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C...	train
1	-0.644026	-0.218808	2.287883	0.424803	-0.721806	-0.398991	-0.414869	0.044067	-0.463779	-0.463790	-0.463790	-0.463751	[H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@...	train
2	-0.459969	0.297628	0.167480	0.130518	0.046509	-0.251763	0.650823	0.217649	-0.113869	-0.113867	-0.113867	-0.113876	[H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]...	train
3	-0.775036	1.089651	-0.104059	1.107464	1.156964	0.315530	1.703451	1.364106	0.237038	0.237068	0.237068	0.237016	[H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H]	test
4	0.346542	-0.820930	-0.902704	0.618991	1.072929	0.839014	-0.128690	0.650888	-1.836564	-1.836554	-1.836554	-1.836579	[H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2(...	train
5	0.409150	-2.145600	-0.758948	-0.494088	-0.113557	-1.271677	-1.391723	-1.357742	-0.503130	-0.503165	-0.503165	-0.503080	[H][N-]C1OC([H])([H])[C@]12[N@@H+]1C([H])([H])...	valid
6	-0.469601	0.294155	0.754484	0.791158	0.412659	-0.396482	0.637064	0.321416	-0.112035	-0.112037	-0.112037	-0.112018	[H]N([H])C1[NH2+][C@]2([H])C([H])([H])[C@]2(C(...	test
7	0.327348	-0.128489	-0.611199	-1.142716	-0.835853	-0.335116	-1.578459	-1.002448	-0.027435	-0.027466	-0.027466	-0.027395	[H]C1([H])[N@@H+]2[C@@]3([H])[C@]2([H])[C@@]2(...	train
8	1.261589	-2.541611	-1.797187	-0.974553	-0.073540	-1.369651	-1.949059	-2.116866	0.891255	0.891221	0.891221	0.891306	[H]C([H])([H])O[C@@]1(C([H])([H])[H])C([H])([H...	train
9	1.653432	-0.118068	0.462979	-1.759313	-1.990326	-0.729647	-1.155817	-1.150208	-0.435730	-0.435769	-0.435769	-0.435686	[H]C(=O)c1c([H])c(N([H])[H])nn1[H]	train

Wizualizacja¶

In [7]:

# ============================================================
# KOMÓRKA 4: WIZUALIZACJA STRUKTUR MOLEKULARNYCH
# ============================================================

from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

print("="*70)
print("WIZUALIZACJA PRZYKŁADOWYCH CZĄSTECZEK")
print("="*70)

# ============================================================
# Wybór przykładowych cząsteczek
# ============================================================

# Losowe próbki z datasetu
n_samples = 12
sample_molecules = df.sample(n=n_samples, random_state=42)

print(f"\n🔬 Wybieram {n_samples} losowych cząsteczek do wizualizacji...\n")

# ============================================================
# Przygotowanie molekuł RDKit
# ============================================================

molecules = []
labels = []

for idx, row in sample_molecules.iterrows():
    smiles = row['SMILES']
    mol = Chem.MolFromSmiles(smiles)

    if mol is not None:
        # Informacje o cząsteczce
        gap = row['gap'] if 'gap' in row else None
        homo = row['HOMO'] if 'HOMO' in row else None

        # Label z kluczowymi informacjami
        label = f"ID: {idx}\n"
        label += f"SMILES: {smiles[:20]}{'...' if len(smiles) > 20 else ''}\n"
        if gap is not None:
            label += f"Gap: {gap:.3f} eV"

        molecules.append(mol)
        labels.append(label)

        print(f"  ✓ Cząsteczka {len(molecules):2d}: {smiles[:40]}{'...' if len(smiles) > 40 else ''}")

print(f"\n✅ Przygotowano {len(molecules)} struktur do wizualizacji")

# ============================================================
# Wizualizacja - Grid 3x4
# ============================================================

print("\n🎨 Generowanie wizualizacji...\n")

fig = plt.figure(figsize=(16, 12))
gs = GridSpec(3, 4, figure=fig, hspace=0.4, wspace=0.3)

for i, (mol, label) in enumerate(zip(molecules, labels)):
    row = i // 4
    col = i % 4

    ax = fig.add_subplot(gs[row, col])

    # Rysowanie cząsteczki
    img = Draw.MolToImage(mol, size=(400, 400))
    ax.imshow(img)
    ax.axis('off')
    ax.set_title(label, fontsize=9, pad=10)

plt.suptitle('Przykładowe cząsteczki z datasetu QM9', 
             fontsize=16, fontweight='bold', y=0.98)

plt.tight_layout()
plt.savefig('qm9_example_molecules.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_example_molecules.png")
plt.show()

# ============================================================
# Wizualizacja cząsteczek o ekstremalnych właściwościach
# ============================================================

print("\n" + "="*70)
print("CZĄSTECZKI O EKSTREMALNYCH WŁAŚCIWOŚCIACH")
print("="*70)

# Najwyższy i najniższy gap
idx_max_gap = df['gap'].idxmax()
idx_min_gap = df['gap'].idxmin()

# Najwyższy i najniższy moment dipolowy
idx_max_dipole = df['dipole_moment'].idxmax()
idx_min_dipole = df['dipole_moment'].idxmin()

extreme_indices = [idx_max_gap, idx_min_gap, idx_max_dipole, idx_min_dipole]
extreme_labels_text = [
    f"Najwyższy gap: {df.loc[idx_max_gap, 'gap']:.3f} eV",
    f"Najniższy gap: {df.loc[idx_min_gap, 'gap']:.3f} eV",
    f"Najwyższy dipole: {df.loc[idx_max_dipole, 'dipole_moment']:.3f} D",
    f"Najniższy dipole: {df.loc[idx_min_dipole, 'dipole_moment']:.3f} D"
]

print("\n🔍 Ekstremalne cząsteczki:")
for i, (idx, text) in enumerate(zip(extreme_indices, extreme_labels_text), 1):
    smiles = df.loc[idx, 'SMILES']
    print(f"  {i}. {text}")
    print(f"     SMILES: {smiles}")

# Wizualizacja ekstremalnych cząsteczek
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
axes = axes.flatten()

for i, (idx, label_text) in enumerate(zip(extreme_indices, extreme_labels_text)):
    smiles = df.loc[idx, 'SMILES']
    mol = Chem.MolFromSmiles(smiles)

    if mol is not None:
        img = Draw.MolToImage(mol, size=(500, 500))
        axes[i].imshow(img)
        axes[i].axis('off')
        axes[i].set_title(f"{label_text}\nSMILES: {smiles[:30]}{'...' if len(smiles) > 30 else ''}", 
                         fontsize=10, fontweight='bold', pad=10)

plt.suptitle('Cząsteczki o ekstremalnych właściwościach', 
             fontsize=14, fontweight='bold', y=0.98)
plt.tight_layout()
plt.savefig('qm9_extreme_molecules.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_extreme_molecules.png")
plt.show()

# ============================================================
# Przykłady według rozmiaru
# ============================================================

print("\n" + "="*70)
print("CZĄSTECZKI WEDŁUG ROZMIARU")
print("="*70)

# Dodanie kolumny z liczbą atomów
df['num_atoms'] = df['SMILES'].apply(lambda x: Chem.MolFromSmiles(x).GetNumAtoms() 
                                      if Chem.MolFromSmiles(x) is not None else 0)

# Najmniejsze i największe
idx_smallest = df['num_atoms'].idxmin()
idx_largest = df['num_atoms'].idxmax()

# Średnie
median_atoms = df['num_atoms'].median()
idx_medium = df.iloc[(df['num_atoms'] - median_atoms).abs().argsort()[:1]].index[0]

size_indices = [idx_smallest, idx_medium, idx_largest]
size_labels = [
    f"Najmniejsza ({df.loc[idx_smallest, 'num_atoms']} atomów)",
    f"Średnia ({df.loc[idx_medium, 'num_atoms']} atomów)",
    f"Największa ({df.loc[idx_largest, 'num_atoms']} atomów)"
]

print("\n📏 Cząsteczki według rozmiaru:")
for label, idx in zip(size_labels, size_indices):
    smiles = df.loc[idx, 'SMILES']
    print(f"  • {label}: {smiles}")

# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for i, (idx, label) in enumerate(zip(size_indices, size_labels)):
    smiles = df.loc[idx, 'SMILES']
    mol = Chem.MolFromSmiles(smiles)

    if mol is not None:
        img = Draw.MolToImage(mol, size=(500, 500))
        axes[i].imshow(img)
        axes[i].axis('off')
        axes[i].set_title(f"{label}\n{smiles}", fontsize=10, fontweight='bold')

plt.suptitle('Porównanie rozmiaru cząsteczek', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_size_comparison.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_size_comparison.png")
plt.show()

print("\n" + "="*70)
print("✅ WIZUALIZACJE ZAKOŃCZONE!")
print("   Wygenerowano 3 pliki PNG z strukturami molekularnymi")
print("="*70)

======================================================================
WIZUALIZACJA PRZYKŁADOWYCH CZĄSTECZEK
======================================================================

🔬 Wybieram 12 losowych cząsteczek do wizualizacji...

  ✓ Cząsteczka  1: [H]Oc1c([H])c(N([H])[H])n([H])c1C([H])([...
  ✓ Cząsteczka  2: [H]O[C@]1([H])C([H])([H])[N@@H+]2C([H])(...
  ✓ Cząsteczka  3: [H]C(=O)[C@]1(C([H])([H])[H])C([H])([H])...
  ✓ Cząsteczka  4: [H]OC([H])([H])[C@](O[H])(C([H])([H])[H]...
  ✓ Cząsteczka  5: [H]OC([H])([H])[C@]12[C@@]3([H])[C@]4([H...
  ✓ Cząsteczka  6: [H][N-]C1O[C@@]2([H])[C@@]([H])(O1)[C@@]...
  ✓ Cząsteczka  7: [H]C([H])([H])OC(=O)[C@@]([H])(OC([H])([...
  ✓ Cząsteczka  8: [H]C1([H])O[C@@]12[C@@]1([H])[N@H+]3[C@]...
  ✓ Cząsteczka  9: [H]OC([H])([H])C([H])([H])N(C([H])=O)C([...
  ✓ Cząsteczka 10: [H]C1=C([H])[C@@]([H])(C([H])([H])[H])OC...
  ✓ Cząsteczka 11: [H]C#C[C@@]1([C@]([H])(O[H])C([H])([H])O...
  ✓ Cząsteczka 12: [H]Oc1nc(N([H])C([H])([H])[H])oc1C([H])(...

✅ Przygotowano 12 struktur do wizualizacji

🎨 Generowanie wizualizacji...

C:\Users\slast\AppData\Local\Temp\ipykernel_6436\3317227650.py:77: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

✅ Zapisano: qm9_example_molecules.png

======================================================================
CZĄSTECZKI O EKSTREMALNYCH WŁAŚCIWOŚCIACH
======================================================================

🔍 Ekstremalne cząsteczki:
  1. Najwyższy gap: 4.752 eV
     SMILES: [H]C#C[C@@]1(C([H])([H])C([H])=O)N([H])[C@]1([H])C([H])([H])[H]
  2. Najniższy gap: -2.537 eV
     SMILES: [H]C1=C([H])c2c(oc([H])c2[H])C1=O
  3. Najwyższy dipole: 4.349 D
     SMILES: [H]C1=C([H])c2c(oc([H])c2[H])C1=O
  4. Najniższy dipole: -1.694 D
     SMILES: [H]C#C[C@@]1(C([H])([H])C([H])=O)N([H])[C@]1([H])C([H])([H])[H]

✅ Zapisano: qm9_extreme_molecules.png

======================================================================
CZĄSTECZKI WEDŁUG ROZMIARU
======================================================================

📏 Cząsteczki według rozmiaru:
  • Najmniejsza (4 atomów): [H]N([H])C(=O)C([H])([H])[H]
  • Średnia (9 atomów): [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H]
  • Największa (11 atomów): [H]/N=C(/C(=O)N(/C([H])=N/[H])C([H])([H])[H])N([H])[H]

✅ Zapisano: qm9_size_comparison.png

======================================================================
✅ WIZUALIZACJE ZAKOŃCZONE!
   Wygenerowano 3 pliki PNG z strukturami molekularnymi
======================================================================

brakujace wartosaci i duplikaty¶

In [8]:

# ============================================================
# SPRAWDZANIE BRAKUJĄCYCH WARTOŚCI
# ============================================================
missing_counts = df.isnull().sum()
total_missing = missing_counts.sum()

if total_missing == 0:
    print("✔️ Brak brakujących wartości w zbiorze danych.")
else:
    print(f"⚠️ Znaleziono {total_missing} brakujących wartości w zbiorze danych:")
    print(missing_counts[missing_counts > 0])

# ============================================================
# SPRAWDZANIE DUPLIKATÓW
# ============================================================
duplicate_count = df.duplicated().sum()

if duplicate_count == 0:
    print("✔️ Brak duplikatów w zbiorze danych.")
else:
    print(f"⚠️ Znaleziono {duplicate_count} duplikatów w zbiorze danych.")

✔️ Brak brakujących wartości w zbiorze danych.
✔️ Brak duplikatów w zbiorze danych.

teraz czas na sprawdzenie czy mamy jkaieś odstające wartosci i wyrysoanie box plotów¶

In [9]:

import matplotlib.pyplot as plt
import seaborn as sns

# ============================================================
# WYBIERAMY TYLKO KOLUMNY NUMERYCZNE
# ============================================================
numeric_cols = df.select_dtypes(include=['float64', 'int64']).columns

print("📊 Sprawdzenie wartości odstających w kolumnach numerycznych:\n")

for col in numeric_cols:
    Q1 = df[col].quantile(0.25)
    Q3 = df[col].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    outliers = df[(df[col] < lower_bound) | (df[col] > upper_bound)]

    print(f"{col}: {len(outliers)} wartości odstających")

# ============================================================
# BOX PLOTY DLA KOLUMN NUMERYCZNYCH
# ============================================================
plt.figure(figsize=(15, len(numeric_cols)*1.5))
for i, col in enumerate(numeric_cols, 1):
    plt.subplot(len(numeric_cols), 1, i)
    sns.boxplot(x=df[col], color='skyblue')
    plt.title(f"Boxplot: {col}")
    plt.tight_layout()

plt.show()

📊 Sprawdzenie wartości odstających w kolumnach numerycznych:

dipole_moment: 10 wartości odstających
polarizability: 33 wartości odstających
HOMO: 68 wartości odstających
LUMO: 0 wartości odstających
gap: 2 wartości odstających
electronic_spatial_extent: 70 wartości odstających
zpve: 3 wartości odstających
U0: 17 wartości odstających
U: 25 wartości odstających
H: 25 wartości odstających
G: 25 wartości odstających
Cv: 25 wartości odstających
num_atoms: 195 wartości odstających

wzory kilku czasteczek o wartosciach odstajacych¶

In [10]:

print(df.columns)

Index(['dipole_moment', 'polarizability', 'HOMO', 'LUMO', 'gap',
       'electronic_spatial_extent', 'zpve', 'U0', 'U', 'H', 'G', 'Cv',
       'SMILES', 'dataset_split', 'num_atoms'],
      dtype='object')

In [11]:

# ============================================================
# WYŚWIETLANIE CZĄSTECZEK Z WARTOŚCIAMI ODSTAJĄCYMI
# ============================================================

from rdkit import Chem
from rdkit.Chem import Draw

# Nazwa kolumny SMILES w Twoim df
smiles_col = 'SMILES'

# Sprawdzenie, które kolumny numeryczne mają wartości odstające
numeric_cols = df.select_dtypes(include=['float64', 'int64']).columns

# Zbieramy wszystkie indeksy z wartościami odstającymi
outlier_indices = set()

for col in numeric_cols:
    Q1 = df[col].quantile(0.25)
    Q3 = df[col].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    outliers = df[(df[col] < lower_bound) | (df[col] > upper_bound)]
    outlier_indices.update(outliers.index.tolist())

# Wybieramy np. pierwsze 5 cząsteczek z wartościami odstającymi
sample_outliers = df.loc[list(outlier_indices)[:5], smiles_col]

# Konwersja SMILES na obiekty RDKit
mols = [Chem.MolFromSmiles(s) for s in sample_outliers]

# Wyświetlenie cząsteczek
Draw.MolsToGridImage(
    mols, 
    molsPerRow=5, 
    subImgSize=(200,200), 
    legends=[f"{i+1}" for i in range(len(mols))]
)

Out[11]:

kolejny punkt¶

teraz czas na korelacja¶

In [12]:

# ============================================================
# KOMÓRKA 5: ANALIZA KORELACJI MIĘDZY WŁAŚCIWOŚCIAMI
# ============================================================

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

print("="*70)
print("ANALIZA KORELACJI - WŁAŚCIWOŚCI QM9")
print("="*70)

# ============================================================
# Przygotowanie danych numerycznych
# ============================================================

# Wybór tylko kolumn numerycznych (bez SMILES i split)
numeric_cols = df.select_dtypes(include=[np.number]).columns.tolist()

# Usunięcie kolumny num_atoms jeśli istnieje (dodana w poprzedniej komórce)
if 'num_atoms' in numeric_cols:
    numeric_cols.remove('num_atoms')

print(f"\n📊 Analizowane właściwości ({len(numeric_cols)}):")
for i, col in enumerate(numeric_cols, 1):
    print(f"  {i:2d}. {col}")

df_numeric = df[numeric_cols].copy()

# ============================================================
# 1. MACIERZ KORELACJI - PEŁNA
# ============================================================

print("\n" + "="*70)
print("1. PEŁNA MACIERZ KORELACJI PEARSONA")
print("="*70)

# Obliczenie korelacji
correlation_matrix = df_numeric.corr(method='pearson')

print(f"\n✓ Macierz korelacji obliczona: {correlation_matrix.shape}")

# Wizualizacja - pełna macierz
plt.figure(figsize=(14, 12))
mask = np.triu(np.ones_like(correlation_matrix, dtype=bool), k=1)

sns.heatmap(correlation_matrix, 
            mask=mask,
            annot=True, 
            fmt='.2f', 
            cmap='RdBu_r', 
            center=0,
            square=True,
            linewidths=0.5,
            cbar_kws={"shrink": 0.8, "label": "Korelacja Pearsona"},
            vmin=-1, 
            vmax=1)

plt.title('Macierz korelacji właściwości kwantowo-chemicznych\n(Korelacja Pearsona)', 
          fontsize=14, fontweight='bold', pad=20)
plt.xticks(rotation=45, ha='right')
plt.yticks(rotation=0)
plt.tight_layout()
plt.savefig('qm9_correlation_full.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_correlation_full.png")
plt.show()

# ============================================================
# 2. NAJSILNIEJSZE KORELACJE
# ============================================================

print("\n" + "="*70)
print("2. TOP 15 NAJSILNIEJSZYCH KORELACJI")
print("="*70)

# Ekstrakcja par korelacji (bez diagonali i duplikatów)
correlations_list = []
for i in range(len(correlation_matrix.columns)):
    for j in range(i+1, len(correlation_matrix.columns)):
        correlations_list.append({
            'Feature 1': correlation_matrix.columns[i],
            'Feature 2': correlation_matrix.columns[j],
            'Correlation': correlation_matrix.iloc[i, j]
        })

df_correlations = pd.DataFrame(correlations_list)
df_correlations['Abs_Correlation'] = df_correlations['Correlation'].abs()
df_correlations = df_correlations.sort_values('Abs_Correlation', ascending=False)

print("\n🔝 TOP 15 najsilniejszych korelacji (wartość bezwzględna):")
print("-" * 70)
for idx, row in df_correlations.head(15).iterrows():
    print(f"  {row['Feature 1']:25s} ↔ {row['Feature 2']:25s}  r = {row['Correlation']:+.4f}")

# Wizualizacja top korelacji
top_n = 15
top_correlations = df_correlations.head(top_n).copy()
top_correlations['Pair'] = (top_correlations['Feature 1'] + '\n↔\n' + 
                             top_correlations['Feature 2'])

plt.figure(figsize=(12, 8))
colors = ['#d73027' if x < 0 else '#4575b4' for x in top_correlations['Correlation']]
bars = plt.barh(range(len(top_correlations)), 
                top_correlations['Correlation'], 
                color=colors, 
                edgecolor='black', 
                linewidth=0.7)

plt.yticks(range(len(top_correlations)), top_correlations['Pair'], fontsize=9)
plt.xlabel('Współczynnik korelacji Pearsona', fontsize=12, fontweight='bold')
plt.title('Top 15 najsilniejszych korelacji między właściwościami', 
          fontsize=14, fontweight='bold', pad=20)
plt.axvline(x=0, color='black', linewidth=0.8, linestyle='-')
plt.grid(axis='x', alpha=0.3, linestyle='--')

# Dodanie wartości na słupkach
for i, (bar, val) in enumerate(zip(bars, top_correlations['Correlation'])):
    plt.text(val + 0.02 if val > 0 else val - 0.02, 
             i, 
             f'{val:.3f}', 
             va='center', 
             ha='left' if val > 0 else 'right',
             fontweight='bold',
             fontsize=9)

plt.tight_layout()
plt.savefig('qm9_top_correlations.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_top_correlations.png")
plt.show()

# ============================================================
# 3. KORELACJA HOMO-LUMO-GAP (kluczowe właściwości)
# ============================================================

print("\n" + "="*70)
print("3. ANALIZA HOMO-LUMO-GAP")
print("="*70)

# Korelacje z gap
gap_correlations = correlation_matrix['gap'].sort_values(ascending=False)
print("\n📊 Korelacje z HOMO-LUMO gap:")
print("-" * 70)
for prop, corr in gap_correlations.items():
    if prop != 'gap':
        print(f"  {prop:30s}  r = {corr:+.4f}")

# Scatter plots - HOMO, LUMO, gap
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Sample dla szybszej wizualizacji
df_sample = df_numeric.sample(n=min(5000, len(df_numeric)), random_state=42)

# HOMO vs LUMO
axes[0].scatter(df_sample['HOMO'], df_sample['LUMO'], 
                alpha=0.3, s=10, c='#1f77b4', edgecolors='none')
axes[0].set_xlabel('HOMO (eV)', fontsize=11, fontweight='bold')
axes[0].set_ylabel('LUMO (eV)', fontsize=11, fontweight='bold')
axes[0].set_title(f'HOMO vs LUMO\nr = {correlation_matrix.loc["HOMO", "LUMO"]:.3f}', 
                  fontsize=12, fontweight='bold')
axes[0].grid(alpha=0.3)

# HOMO vs gap
axes[1].scatter(df_sample['HOMO'], df_sample['gap'], 
                alpha=0.3, s=10, c='#ff7f0e', edgecolors='none')
axes[1].set_xlabel('HOMO (eV)', fontsize=11, fontweight='bold')
axes[1].set_ylabel('gap (eV)', fontsize=11, fontweight='bold')
axes[1].set_title(f'HOMO vs gap\nr = {correlation_matrix.loc["HOMO", "gap"]:.3f}', 
                  fontsize=12, fontweight='bold')
axes[1].grid(alpha=0.3)

# LUMO vs gap
axes[2].scatter(df_sample['LUMO'], df_sample['gap'], 
                alpha=0.3, s=10, c='#2ca02c', edgecolors='none')
axes[2].set_xlabel('LUMO (eV)', fontsize=11, fontweight='bold')
axes[2].set_ylabel('gap (eV)', fontsize=11, fontweight='bold')
axes[2].set_title(f'LUMO vs gap\nr = {correlation_matrix.loc["LUMO", "gap"]:.3f}', 
                  fontsize=12, fontweight='bold')
axes[2].grid(alpha=0.3)

plt.suptitle('Zależności między kluczowymi właściwościami elektronowymi', 
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('qm9_homo_lumo_gap.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_homo_lumo_gap.png")
plt.show()

# ============================================================
# 4. PAIRPLOT - WYBRANE WŁAŚCIWOŚCI
# ============================================================

print("\n" + "="*70)
print("4. PAIRPLOT - KLUCZOWE WŁAŚCIWOŚCI")
print("="*70)

# Wybór najważniejszych właściwości
key_properties = ['HOMO', 'LUMO', 'gap', 'dipole_moment', 'polarizability']

print(f"\n📈 Tworzenie pairplot dla {len(key_properties)} właściwości...")
print(f"   Używam {min(2000, len(df))} próbek dla szybszej wizualizacji")

df_pairplot = df[key_properties].sample(n=min(2000, len(df)), random_state=42)

# Pairplot
pairplot_fig = sns.pairplot(df_pairplot, 
                             diag_kind='kde',
                             plot_kws={'alpha': 0.4, 's': 15, 'edgecolor': 'none'},
                             diag_kws={'linewidth': 2})

pairplot_fig.fig.suptitle('Pairplot kluczowych właściwości kwantowo-chemicznych', 
                          fontsize=14, fontweight='bold', y=1.01)

plt.savefig('qm9_pairplot.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_pairplot.png")
plt.show()

# ============================================================
# 5. KORELACJA SPEARMANA (dla sprawdzenia nieliniowych zależności)
# ============================================================

print("\n" + "="*70)
print("5. PORÓWNANIE KORELACJI PEARSONA vs SPEARMANA")
print("="*70)

# Korelacja Spearmana
spearman_matrix = df_numeric.corr(method='spearman')

# Różnica między Pearson a Spearman
diff_matrix = np.abs(correlation_matrix - spearman_matrix)

print("\n📊 Największe różnice między Pearson a Spearman (nieliniowości):")
print("-" * 70)

# Znalezienie największych różnic
max_diffs = []
for i in range(len(diff_matrix.columns)):
    for j in range(i+1, len(diff_matrix.columns)):
        max_diffs.append({
            'Feature 1': diff_matrix.columns[i],
            'Feature 2': diff_matrix.columns[j],
            'Pearson': correlation_matrix.iloc[i, j],
            'Spearman': spearman_matrix.iloc[i, j],
            'Difference': diff_matrix.iloc[i, j]
        })

df_diffs = pd.DataFrame(max_diffs).sort_values('Difference', ascending=False)

for idx, row in df_diffs.head(10).iterrows():
    print(f"  {row['Feature 1']:20s} ↔ {row['Feature 2']:20s}")
    print(f"    Pearson: {row['Pearson']:+.4f}  |  Spearman: {row['Spearman']:+.4f}  |  Diff: {row['Difference']:.4f}")

# ============================================================
# 6. PODSUMOWANIE STATYSTYCZNE
# ============================================================

print("\n" + "="*70)
print("6. PODSUMOWANIE ANALIZY KORELACJI")
print("="*70)

print(f"""
📊 STATYSTYKI KORELACJI:

• Liczba analizowanych właściwości: {len(numeric_cols)}
• Liczba par korelacji: {len(correlations_list)}

• Najsilniejsza korelacja dodatnia:
  {df_correlations.iloc[0]['Feature 1']} ↔ {df_correlations.iloc[0]['Feature 2']}
  r = {df_correlations.iloc[0]['Correlation']:.4f}

• Najsilniejsza korelacja ujemna:
  {df_correlations[df_correlations['Correlation'] < 0].iloc[0]['Feature 1']} ↔ {df_correlations[df_correlations['Correlation'] < 0].iloc[0]['Feature 2']}
  r = {df_correlations[df_correlations['Correlation'] < 0].iloc[0]['Correlation']:.4f}

• Silne korelacje (|r| > 0.7): {len(df_correlations[df_correlations['Abs_Correlation'] > 0.7])}
• Średnie korelacje (0.3 < |r| < 0.7): {len(df_correlations[(df_correlations['Abs_Correlation'] > 0.3) & (df_correlations['Abs_Correlation'] <= 0.7)])}
• Słabe korelacje (|r| < 0.3): {len(df_correlations[df_correlations['Abs_Correlation'] <= 0.3])}
""")

print("="*70)
print("✅ ANALIZA KORELACJI ZAKOŃCZONA!")
print("   Wygenerowano 4 pliki PNG z wizualizacjami")
print("="*70)

======================================================================
ANALIZA KORELACJI - WŁAŚCIWOŚCI QM9
======================================================================

📊 Analizowane właściwości (12):
   1. dipole_moment
   2. polarizability
   3. HOMO
   4. LUMO
   5. gap
   6. electronic_spatial_extent
   7. zpve
   8. U0
   9. U
  10. H
  11. G
  12. Cv

======================================================================
1. PEŁNA MACIERZ KORELACJI PEARSONA
======================================================================

✓ Macierz korelacji obliczona: (12, 12)
✅ Zapisano: qm9_correlation_full.png

======================================================================
2. TOP 15 NAJSILNIEJSZYCH KORELACJI
======================================================================

🔝 TOP 15 najsilniejszych korelacji (wartość bezwzględna):
----------------------------------------------------------------------
  H                         ↔ G                          r = +1.0000
  U                         ↔ H                          r = +1.0000
  U                         ↔ G                          r = +1.0000
  U                         ↔ Cv                         r = +1.0000
  H                         ↔ Cv                         r = +1.0000
  G                         ↔ Cv                         r = +1.0000
  LUMO                      ↔ gap                        r = +0.8670
  electronic_spatial_extent ↔ U0                         r = +0.8288
  polarizability            ↔ U0                         r = +0.7982
  zpve                      ↔ U0                         r = +0.7758
  polarizability            ↔ zpve                       r = +0.7655
  polarizability            ↔ electronic_spatial_extent  r = +0.7170
  LUMO                      ↔ zpve                       r = +0.6705
  gap                       ↔ zpve                       r = +0.5705
  electronic_spatial_extent ↔ zpve                       r = +0.5284

✅ Zapisano: qm9_top_correlations.png

======================================================================
3. ANALIZA HOMO-LUMO-GAP
======================================================================

📊 Korelacje z HOMO-LUMO gap:
----------------------------------------------------------------------
  LUMO                            r = +0.8670
  zpve                            r = +0.5705
  G                               r = +0.3873
  H                               r = +0.3873
  U                               r = +0.3873
  Cv                              r = +0.3873
  U0                              r = +0.2404
  polarizability                  r = +0.1883
  electronic_spatial_extent       r = +0.0438
  HOMO                            r = -0.2878
  dipole_moment                   r = -0.3407

✅ Zapisano: qm9_homo_lumo_gap.png

======================================================================
4. PAIRPLOT - KLUCZOWE WŁAŚCIWOŚCI
======================================================================

📈 Tworzenie pairplot dla 5 właściwości...
   Używam 1000 próbek dla szybszej wizualizacji
✅ Zapisano: qm9_pairplot.png

======================================================================
5. PORÓWNANIE KORELACJI PEARSONA vs SPEARMANA
======================================================================

📊 Największe różnice między Pearson a Spearman (nieliniowości):
----------------------------------------------------------------------
  polarizability       ↔ H                   
    Pearson: -0.1031  |  Spearman: +0.1497  |  Diff: 0.2528
  polarizability       ↔ G                   
    Pearson: -0.1031  |  Spearman: +0.1497  |  Diff: 0.2528
  polarizability       ↔ U                   
    Pearson: -0.1031  |  Spearman: +0.1496  |  Diff: 0.2527
  polarizability       ↔ Cv                  
    Pearson: -0.1031  |  Spearman: +0.1495  |  Diff: 0.2526
  U0                   ↔ H                   
    Pearson: -0.3127  |  Spearman: -0.1948  |  Diff: 0.1179
  U0                   ↔ G                   
    Pearson: -0.3127  |  Spearman: -0.1948  |  Diff: 0.1179
  U0                   ↔ U                   
    Pearson: -0.3127  |  Spearman: -0.1950  |  Diff: 0.1177
  U0                   ↔ Cv                  
    Pearson: -0.3128  |  Spearman: -0.1953  |  Diff: 0.1175
  polarizability       ↔ gap                 
    Pearson: +0.1883  |  Spearman: +0.2940  |  Diff: 0.1057
  HOMO                 ↔ U0                  
    Pearson: +0.1213  |  Spearman: +0.0184  |  Diff: 0.1028

======================================================================
6. PODSUMOWANIE ANALIZY KORELACJI
======================================================================

📊 STATYSTYKI KORELACJI:

• Liczba analizowanych właściwości: 12
• Liczba par korelacji: 66

• Najsilniejsza korelacja dodatnia:
  H ↔ G
  r = 1.0000

• Najsilniejsza korelacja ujemna:
  dipole_moment ↔ LUMO
  r = -0.3941

• Silne korelacje (|r| > 0.7): 12
• Średnie korelacje (0.3 < |r| < 0.7): 24
• Słabe korelacje (|r| < 0.3): 30

======================================================================
✅ ANALIZA KORELACJI ZAKOŃCZONA!
   Wygenerowano 4 pliki PNG z wizualizacjami
======================================================================

🔝 TOP 15 NAJSILNIEJSZYCH KORELACJI¶

Poniżej przedstawiono 15 par właściwości z najsilniejszymi korelacjami w datasecie QM9:

Korelacje perfekcyjne (r = 1.00)¶

Te właściwości są ze sobą matematycznie związane poprzez definicje termodynamiczne:

• H ↔ G → r = +1.0000 (Entalpia vs Energia Gibbsa)
• U ↔ H → r = +1.0000 (Energia wewnętrzna vs Entalpia)
• U ↔ G → r = +1.0000 (Energia wewnętrzna vs Energia Gibbsa)
• U ↔ Cv → r = +1.0000 (Energia wewnętrzna vs Pojemność cieplna)
• H ↔ Cv → r = +1.0000 (Entalpia vs Pojemność cieplna)
• G ↔ Cv → r = +1.0000 (Energia Gibbsa vs Pojemność cieplna)

⚠️ Uwaga: Te 6 par wykazują korelację 1.0, co sugeruje, że są ze sobą liniowo zależne. W modelowaniu ML należy rozważyć usunięcie redundantnych zmiennych, aby uniknąć multikolinearności.

Korelacje bardzo silne (r > 0.70)¶

• LUMO ↔ gap → r = +0.8670
  LUMO energy jest głównym składnikiem HOMO-LUMO gap

• electronic_spatial_extent ↔ U0 → r = +0.8288
  Większe cząsteczki mają wyższą energię wewnętrzną

• polarizability ↔ U0 → r = +0.7982
  Polaryzowalność rośnie z rozmiarem molekuły

• zpve ↔ U0 → r = +0.7758
  Energia drgań punktu zerowego koreluje z całkowitą energią

• polarizability ↔ zpve → r = +0.7655
  Obie właściwości zależą od rozmiaru molekuły

• polarizability ↔ electronic_spatial_extent → r = +0.7170
  Bezpośrednia zależność od objętości elektronowej

Korelacje silne (0.50 < r < 0.70)¶

• LUMO ↔ zpve → r = +0.6705
• gap ↔ zpve → r = +0.5705
• electronic_spatial_extent ↔ zpve → r = +0.5284

---

💡 Kluczowe wnioski:¶

1. Właściwości termodynamiczne (U, H, G, Cv) są ze sobą idealnie skorelowane - wystarczy jedna z nich w modelach predykcyjnych
2. LUMO i gap są silnie powiązane - gap częściowo zależy od LUMO
3. Rozmiar molekuły (mierzony przez polarizability i electronic_spatial_extent) silnie wpływa na właściwości energetyczne
4. Brak silnych korelacji ujemnych w top 15 - większość zależności jest dodatnia

Walidacja SMILES - ile błędnych struktur, dlaczego¶

In [13]:

# ============================================================
# KOMÓRKA 7: WALIDACJA I ANALIZA STRUKTUR SMILES
# ============================================================

from rdkit import Chem, RDLogger
from rdkit.Chem import Descriptors, rdMolDescriptors
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter

# Wyciszenie ostrzeżeń RDKit
RDLogger.DisableLog('rdApp.*')

print("="*70)
print("WALIDACJA STRUKTUR SMILES - QM9 DATASET")
print("="*70)

# ============================================================
# 1. PODSTAWOWA WALIDACJA SMILES
# ============================================================

print("\n📋 ETAP 1: PODSTAWOWA WALIDACJA")
print("-" * 70)

def validate_smiles_detailed(smiles):
    """
    Szczegółowa walidacja SMILES z diagnozą błędów
    
    Returns:
        dict: {'valid': bool, 'error': str, 'mol': Mol object}
    """
    try:
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            return {'valid': False, 'error': 'Parse Error', 'mol': None}

        # Próba sanityzacji
        Chem.SanitizeMol(mol)
        return {'valid': True, 'error': None, 'mol': mol}

    except Exception as e:
        error_type = str(e).split(':')[0] if ':' in str(e) else str(e)
        return {'valid': False, 'error': error_type, 'mol': None}

# Walidacja wszystkich SMILES
print("Walidacja struktur SMILES...")
validation_results = df['SMILES'].apply(validate_smiles_detailed)

df['is_valid'] = validation_results.apply(lambda x: x['valid'])
df['error_type'] = validation_results.apply(lambda x: x['error'])

# Statystyki
total_smiles = len(df)
valid_smiles = df['is_valid'].sum()
invalid_smiles = total_smiles - valid_smiles
validity_pct = (valid_smiles / total_smiles) * 100

print(f"\n✅ Poprawne SMILES: {valid_smiles:,} ({validity_pct:.2f}%)")
print(f"❌ Niepoprawne SMILES: {invalid_smiles:,} ({(100-validity_pct):.2f}%)")

# ============================================================
# 2. ANALIZA BŁĘDÓW WALIDACJI
# ============================================================

if invalid_smiles > 0:
    print("\n" + "="*70)
    print("📊 ANALIZA BŁĘDÓW WALIDACJI")
    print("="*70)

    # Typy błędów
    error_counts = df[~df['is_valid']]['error_type'].value_counts()

    print("\n🔍 Typy błędów walidacji:")
    for error, count in error_counts.items():
        pct = (count / invalid_smiles) * 100
        print(f"  • {error:40s}  {count:5d} ({pct:5.2f}%)")

    # Przykłady błędnych SMILES
    print("\n⚠️  Przykłady niepoprawnych SMILES:")
    invalid_samples = df[~df['is_valid']].head(10)
    for idx, row in invalid_samples.iterrows():
        print(f"\n  [{idx}] {row['SMILES'][:60]}{'...' if len(row['SMILES']) > 60 else ''}")
        print(f"      Błąd: {row['error_type']}")

    # Wizualizacja typów błędów
    if len(error_counts) > 0:
        plt.figure(figsize=(12, 6))
        colors = plt.cm.Reds(np.linspace(0.4, 0.8, len(error_counts)))
        plt.barh(range(len(error_counts)), error_counts.values, color=colors, edgecolor='black')
        plt.yticks(range(len(error_counts)), error_counts.index)
        plt.xlabel('Liczba wystąpień', fontsize=12, fontweight='bold')
        plt.title('Typy błędów walidacji SMILES', fontsize=14, fontweight='bold')
        plt.grid(axis='x', alpha=0.3)

        for i, val in enumerate(error_counts.values):
            plt.text(val + max(error_counts.values)*0.01, i, f'{val}', 
                    va='center', fontweight='bold')

        plt.tight_layout()
        plt.savefig('qm9_smiles_errors.png', dpi=300, bbox_inches='tight')
        print("\n✅ Zapisano: qm9_smiles_errors.png")
        plt.show()

else:
    print("\n✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.")

# ============================================================
# 3. ANALIZA SKŁADU CHEMICZNEGO
# ============================================================

print("\n" + "="*70)
print("🧪 ANALIZA SKŁADU CHEMICZNEGO")
print("="*70)

# Tylko dla poprawnych SMILES
df_valid = df[df['is_valid']].copy()

def analyze_composition(smiles):
    """Analiza składu atomowego cząsteczki"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    atoms = [atom.GetSymbol() for atom in mol.GetAtoms()]
    return Counter(atoms)

print("\nAnalizuję skład atomowy cząsteczek...")
compositions = df_valid['SMILES'].apply(analyze_composition)

# Ekstrakcja wszystkich typów atomów
all_atoms = []
for comp in compositions.dropna():
    all_atoms.extend(comp.keys())

atom_types = Counter(all_atoms)

print(f"\n📊 Typy atomów występujące w datasecie ({len(atom_types)}):")
for atom, count in atom_types.most_common():
    pct = (count / len(df_valid)) * 100
    print(f"  • {atom:3s}  występuje w {count:6,} cząsteczkach ({pct:5.2f}%)")

# Rozkład liczby atomów per typ
atom_counts = {atom: [] for atom in atom_types.keys()}
for comp in compositions.dropna():
    for atom in atom_types.keys():
        atom_counts[atom].append(comp.get(atom, 0))

# Konwersja do DataFrame
df_atoms = pd.DataFrame(atom_counts)

print(f"\n📈 Statystyki liczby atomów:")
print(df_atoms.describe())

# ============================================================
# 4. WIZUALIZACJA SKŁADU ATOMOWEGO
# ============================================================

print("\n" + "="*70)
print("📊 WIZUALIZACJA SKŁADU ATOMOWEGO")
print("="*70)

# Wykresy dla głównych atomów (H, C, N, O, F)
main_atoms = ['H', 'C', 'N', 'O', 'F']
main_atoms = [atom for atom in main_atoms if atom in df_atoms.columns]

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()

for i, atom in enumerate(main_atoms):
    if atom in df_atoms.columns:
        data = df_atoms[atom]
        axes[i].hist(data, bins=50, color='steelblue', edgecolor='black', alpha=0.7)
        axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
        axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
        axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.1f}, Max: {data.max():.0f}', 
                         fontsize=12, fontweight='bold')
        axes[i].grid(axis='y', alpha=0.3)

# Wykres porównawczy - średnie liczby atomów
if len(main_atoms) < 6:
    ax_comparison = axes[5]
    means = [df_atoms[atom].mean() for atom in main_atoms]
    colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'][:len(main_atoms)]
    bars = ax_comparison.bar(main_atoms, means, color=colors, edgecolor='black', linewidth=1.5)
    ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
    ax_comparison.set_title('Średni skład atomowy cząsteczek QM9', fontsize=12, fontweight='bold')
    ax_comparison.grid(axis='y', alpha=0.3)

    for bar, val in zip(bars, means):
        ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.1, 
                          f'{val:.1f}', ha='center', fontweight='bold')

plt.tight_layout()
plt.savefig('qm9_atomic_composition.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_atomic_composition.png")
plt.show()

# ============================================================
# 5. ANALIZA WŁAŚCIWOŚCI MOLEKULARNYCH (RDKit)
# ============================================================

print("\n" + "="*70)
print("⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH")
print("="*70)

def calculate_rdkit_descriptors(smiles):
    """Oblicza podstawowe deskryptory RDKit"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    try:
        return {
            'MolWt': Descriptors.MolWt(mol),
            'NumAtoms': mol.GetNumAtoms(),
            'NumHeavyAtoms': mol.GetNumHeavyAtoms(),
            'NumHeteroatoms': rdMolDescriptors.CalcNumHeteroatoms(mol),
            'NumRotatableBonds': rdMolDescriptors.CalcNumRotatableBonds(mol),
            'NumHBD': rdMolDescriptors.CalcNumHBD(mol),
            'NumHBA': rdMolDescriptors.CalcNumHBA(mol),
            'NumRings': rdMolDescriptors.CalcNumRings(mol),
            'NumAromaticRings': rdMolDescriptors.CalcNumAromaticRings(mol),
            'TPSA': Descriptors.TPSA(mol),
            'LogP': Descriptors.MolLogP(mol),
            'FractionCSP3': rdMolDescriptors.CalcFractionCsp3(mol)
        }
    except:
        return None

print("Obliczam deskryptory RDKit dla poprawnych struktur...")
descriptors = df_valid['SMILES'].apply(calculate_rdkit_descriptors)
df_descriptors = descriptors.apply(pd.Series)

# Sprawdzenie czy są jakiekolwiek dane
if df_descriptors.empty or df_descriptors.shape[1] == 0:
    print("\n⚠️  BRAK POPRAWNYCH STRUKTUR DO ANALIZY!")
    print("   Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów.")
else:
    # Dodanie do głównego DataFrame
    df_valid = pd.concat([df_valid.reset_index(drop=True), 
                          df_descriptors.reset_index(drop=True)], axis=1)

    print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur")

    print("\n📊 Statystyki deskryptorów molekularnych:")
    print(df_descriptors.describe())

# ============================================================
# 6. CHARAKTERYSTYKA DATASETU QM9
# ============================================================

if not df_descriptors.empty and df_descriptors.shape[1] > 0:
    print("\n" + "="*70)
    print("📋 CHARAKTERYSTYKA DATASETU QM9")
    print("="*70)

    print(f"""
🔬 SKŁAD CHEMICZNY:
  • Dozwolone atomy: {', '.join(sorted(atom_types.keys()))}
  • Najczęstsze atomy: H, C, N, O, F (QM9 spec)
  • Średnia liczba atomów: {df_descriptors['NumAtoms'].mean():.1f}
  • Średnia liczba ciężkich atomów: {df_descriptors['NumHeavyAtoms'].mean():.1f}
  • Max ciężkich atomów: {df_descriptors['NumHeavyAtoms'].max():.0f} (limit QM9: 9)

📏 ROZMIAR MOLEKUŁ:
  • Masa molowa: {df_descriptors['MolWt'].mean():.1f} ± {df_descriptors['MolWt'].std():.1f} g/mol
  • Zakres: {df_descriptors['MolWt'].min():.1f} - {df_descriptors['MolWt'].max():.1f} g/mol

🔗 STRUKTURY:
  • Średnia liczba pierścieni: {df_descriptors['NumRings'].mean():.2f}
  • Cząsteczki aromatyczne: {(df_descriptors['NumAromaticRings'] > 0).sum():,} ({(df_descriptors['NumAromaticRings'] > 0).sum()/len(df_descriptors)*100:.1f}%)
  • Średnia rotowalne wiązania: {df_descriptors['NumRotatableBonds'].mean():.2f}

💧 WŁAŚCIWOŚCI FIZYKOCHEMICZNE:
  • Średnie LogP: {df_descriptors['LogP'].mean():.2f} ± {df_descriptors['LogP'].std():.2f}
  • Średnie TPSA: {df_descriptors['TPSA'].mean():.1f} Ų
  • Średnia frakcja sp3: {df_descriptors['FractionCSP3'].mean():.2f}

✅ JAKOŚĆ DANYCH:
  • Validacja SMILES: {validity_pct:.2f}%
  • Kompletność: {(df['is_valid'].sum() / len(df) * 100):.2f}%
  • Błędne struktury: {invalid_smiles:,}
""")

# ============================================================
# 7. WIZUALIZACJA PORÓWNAWCZA
# ============================================================

if not df_descriptors.empty and df_descriptors.shape[1] > 0:
    print("\n" + "="*70)
    print("📊 WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH")
    print("="*70)

    fig, axes = plt.subplots(2, 3, figsize=(18, 10))

    # 1. Masa molowa
    axes[0, 0].hist(df_descriptors['MolWt'], bins=50, color='#3498db', 
                    edgecolor='black', alpha=0.7)
    axes[0, 0].set_xlabel('Masa molowa (g/mol)', fontweight='bold')
    axes[0, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 0].set_title('Rozkład masy molowej', fontweight='bold', fontsize=12)
    axes[0, 0].grid(axis='y', alpha=0.3)

    # 2. Liczba ciężkich atomów
    axes[0, 1].hist(df_descriptors['NumHeavyAtoms'], bins=range(1, 11), 
                    color='#e74c3c', edgecolor='black', alpha=0.7)
    axes[0, 1].set_xlabel('Liczba ciężkich atomów', fontweight='bold')
    axes[0, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 1].set_title('Rozkład liczby ciężkich atomów (max=9)', fontweight='bold', fontsize=12)
    axes[0, 1].grid(axis='y', alpha=0.3)

    # 3. LogP
    axes[0, 2].hist(df_descriptors['LogP'], bins=50, color='#2ecc71', 
                    edgecolor='black', alpha=0.7)
    axes[0, 2].set_xlabel('LogP', fontweight='bold')
    axes[0, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 2].set_title('Rozkład lipofilowości (LogP)', fontweight='bold', fontsize=12)
    axes[0, 2].grid(axis='y', alpha=0.3)

    # 4. TPSA
    axes[1, 0].hist(df_descriptors['TPSA'], bins=50, color='#f39c12', 
                    edgecolor='black', alpha=0.7)
    axes[1, 0].set_xlabel('TPSA (Ų)', fontweight='bold')
    axes[1, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 0].set_title('Rozkład polarnej powierzchni (TPSA)', fontweight='bold', fontsize=12)
    axes[1, 0].grid(axis='y', alpha=0.3)

    # 5. Liczba pierścieni
    ring_counts = df_descriptors['NumRings'].value_counts().sort_index()
    axes[1, 1].bar(ring_counts.index, ring_counts.values, color='#9b59b6', 
                   edgecolor='black', alpha=0.7)
    axes[1, 1].set_xlabel('Liczba pierścieni', fontweight='bold')
    axes[1, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 1].set_title('Rozkład liczby pierścieni', fontweight='bold', fontsize=12)
    axes[1, 1].grid(axis='y', alpha=0.3)

    # 6. Frakcja sp3
    axes[1, 2].hist(df_descriptors['FractionCSP3'], bins=50, color='#1abc9c', 
                    edgecolor='black', alpha=0.7)
    axes[1, 2].set_xlabel('Frakcja węgli sp³', fontweight='bold')
    axes[1, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 2].set_title('Rozkład frakcji sp³ (saturacja)', fontweight='bold', fontsize=12)
    axes[1, 2].grid(axis='y', alpha=0.3)

    plt.suptitle('Właściwości molekularne datasetu QM9', fontsize=16, fontweight='bold', y=1.00)
    plt.tight_layout()
    plt.savefig('qm9_molecular_properties.png', dpi=300, bbox_inches='tight')
    print("✅ Zapisano: qm9_molecular_properties.png")
    plt.show()
else:
    print("\n⚠️  Brak danych do wizualizacji właściwości molekularnych")

print("\n" + "="*70)
print("✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!")
print("="*70)

======================================================================
WALIDACJA STRUKTUR SMILES - QM9 DATASET
======================================================================

📋 ETAP 1: PODSTAWOWA WALIDACJA
----------------------------------------------------------------------
Walidacja struktur SMILES...

✅ Poprawne SMILES: 1,000 (100.00%)
❌ Niepoprawne SMILES: 0 (0.00%)

✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.

======================================================================
🧪 ANALIZA SKŁADU CHEMICZNEGO
======================================================================

Analizuję skład atomowy cząsteczek...

📊 Typy atomów występujące w datasecie (5):
  • C    występuje w  1,000 cząsteczkach (100.00%)
  • O    występuje w    843 cząsteczkach (84.30%)
  • N    występuje w    622 cząsteczkach (62.20%)
  • H    występuje w     16 cząsteczkach ( 1.60%)
  • F    występuje w     14 cząsteczkach ( 1.40%)

📈 Statystyki liczby atomów:
                 C            N            O            H            F
count  1000.000000  1000.000000  1000.000000  1000.000000  1000.000000
mean      6.251000     1.093000     1.392000     0.017000     0.020000
std       1.282046     1.138265     0.875844     0.136857     0.188774
min       2.000000     0.000000     0.000000     0.000000     0.000000
25%       5.000000     0.000000     1.000000     0.000000     0.000000
50%       6.000000     1.000000     1.000000     0.000000     0.000000
75%       7.000000     2.000000     2.000000     0.000000     0.000000
max       9.000000     6.000000     4.000000     2.000000     3.000000

======================================================================
📊 WIZUALIZACJA SKŁADU ATOMOWEGO
======================================================================
✅ Zapisano: qm9_atomic_composition.png

======================================================================
⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH
======================================================================
Obliczam deskryptory RDKit dla poprawnych struktur...

⚠️  BRAK POPRAWNYCH STRUKTUR DO ANALIZY!
   Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów.

⚠️  Brak danych do wizualizacji właściwości molekularnych

======================================================================
✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!
======================================================================

In [14]:

# debudding

In [15]:

# ============================================================
# DEBUGGING - SPRAWDZENIE KOLUMNY SMILES
# ============================================================

print("="*70)
print("🔍 DEBUGGING - ANALIZA KOLUMNY SMILES")
print("="*70)

# 1. Sprawdź nazwy kolumn
print("\n1️⃣ Kolumny w DataFrame:")
print(df.columns.tolist())

# 2. Sprawdź czy kolumna SMILES istnieje
print(f"\n2️⃣ Czy 'SMILES' w kolumnach? {('SMILES' in df.columns)}")

# 3. Jeśli istnieje, pokaż przykłady
if 'SMILES' in df.columns:
    print("\n3️⃣ Pierwsze 10 wartości w kolumnie SMILES:")
    for i, smiles in enumerate(df['SMILES'].head(10), 1):
        print(f"  {i:2d}. {smiles}")

    # 4. Typ danych
    print(f"\n4️⃣ Typ danych kolumny SMILES: {df['SMILES'].dtype}")

    # 5. Długość SMILES
    print(f"\n5️⃣ Statystyki długości SMILES:")
    smiles_lengths = df['SMILES'].astype(str).str.len()
    print(f"  Min: {smiles_lengths.min()}")
    print(f"  Max: {smiles_lengths.max()}")
    print(f"  Średnia: {smiles_lengths.mean():.1f}")

    # 6. Czy są puste wartości?
    print(f"\n6️⃣ Puste/NaN wartości: {df['SMILES'].isna().sum()}")

    # 7. Test walidacji pojedynczego SMILES
    print("\n7️⃣ Test walidacji pojedynczego SMILES:")
    test_smiles = df['SMILES'].iloc[0]
    print(f"  SMILES: {test_smiles}")
    print(f"  Typ: {type(test_smiles)}")

    from rdkit import Chem
    mol = Chem.MolFromSmiles(str(test_smiles))
    if mol is None:
        print(f"  ❌ RDKit NIE MOŻE SPARSOWAĆ!")
    else:
        print(f"  ✅ RDKit sparsował poprawnie")
        print(f"  Liczba atomów: {mol.GetNumAtoms()}")

else:
    print("\n❌ KOLUMNA 'SMILES' NIE ISTNIEJE!")
    print("\nMożliwe alternatywne nazwy:")
    for col in df.columns:
        if 'smile' in col.lower() or 'mol' in col.lower() or 'id' in col.lower():
            print(f"  • {col}")

print("\n" + "="*70)

======================================================================
🔍 DEBUGGING - ANALIZA KOLUMNY SMILES
======================================================================

1️⃣ Kolumny w DataFrame:
['dipole_moment', 'polarizability', 'HOMO', 'LUMO', 'gap', 'electronic_spatial_extent', 'zpve', 'U0', 'U', 'H', 'G', 'Cv', 'SMILES', 'dataset_split', 'num_atoms', 'is_valid', 'error_type']

2️⃣ Czy 'SMILES' w kolumnach? True

3️⃣ Pierwsze 10 wartości w kolumnie SMILES:
   1. [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H]
   2. [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H]
   3. [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]1([H])OC([H])([H])[C@@]12C([H])([H])[H]
   4. [H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H]
   5. [H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2([H])[H]
   6. [H][N-]C1OC([H])([H])[C@]12[N@@H+]1C([H])([H])[C@]1([H])C2([H])[H]
   7. [H]N([H])C1[NH2+][C@]2([H])C([H])([H])[C@]2(C([H])([H])C([H])([H])[H])O1
   8. [H]C1([H])[N@@H+]2[C@@]3([H])[C@]2([H])[C@@]2([H])C([H])([H])[C@]1([H])[C@@]32[H]
   9. [H]C([H])([H])O[C@@]1(C([H])([H])[H])C([H])([H])[C@@]2([H])[C@@]([H])(C([H])([H])[H])[C@@]12[H]
  10. [H]C(=O)c1c([H])c(N([H])[H])nn1[H]

4️⃣ Typ danych kolumny SMILES: object

5️⃣ Statystyki długości SMILES:
  Min: 13
  Max: 101
  Średnia: 62.8

6️⃣ Puste/NaN wartości: 0

7️⃣ Test walidacji pojedynczego SMILES:
  SMILES: [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H]
  Typ: <class 'str'>
  ✅ RDKit sparsował poprawnie
  Liczba atomów: 9

======================================================================

SMILES są poprawne, ale mają jawne atomy wodoru [H] i sterochemię [C@], [C@@] - to są bardzo szczegółowe SMILES z QM9! Problem jest w poprzednim kodzie - prawdopodobnie kolumny is_valid i error_type już istnieją z poprzedniego uruchomienia i nadpisują nowe wyniki. Rozwiązanie: Usuńmy te kolumny przed walidacją:

In [16]:

# ============================================================
# KOMÓRKA 7: WALIDACJA I ANALIZA STRUKTUR SMILES
# ============================================================

from rdkit import Chem, RDLogger
from rdkit.Chem import Descriptors, rdMolDescriptors
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter

# Wyciszenie ostrzeżeń RDKit
RDLogger.DisableLog('rdApp.*')

print("="*70)
print("WALIDACJA STRUKTUR SMILES - QM9 DATASET")
print("="*70)

# ============================================================
# 1. PODSTAWOWA WALIDACJA SMILES
# ============================================================

print("\n📋 ETAP 1: PODSTAWOWA WALIDACJA")
print("-" * 70)

# Usunięcie starych kolumn walidacji jeśli istnieją
if 'is_valid' in df.columns:
    df = df.drop(columns=['is_valid'])
if 'error_type' in df.columns:
    df = df.drop(columns=['error_type'])

def validate_smiles_detailed(smiles):
    """
    Szczegółowa walidacja SMILES z diagnozą błędów
    
    Returns:
        dict: {'valid': bool, 'error': str, 'mol': Mol object}
    """
    try:
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            return {'valid': False, 'error': 'Parse Error', 'mol': None}

        # Próba sanityzacji
        Chem.SanitizeMol(mol)
        return {'valid': True, 'error': None, 'mol': mol}

    except Exception as e:
        error_type = str(e).split(':')[0] if ':' in str(e) else str(e)
        return {'valid': False, 'error': error_type, 'mol': None}

# Walidacja wszystkich SMILES
print("Walidacja struktur SMILES...")
print(f"Przykładowe SMILES do walidacji:")
print(df['SMILES'].head(5).tolist())

validation_results = df['SMILES'].apply(validate_smiles_detailed)

df['is_valid'] = validation_results.apply(lambda x: x['valid'])
df['error_type'] = validation_results.apply(lambda x: x['error'])

# Statystyki
total_smiles = len(df)
valid_smiles = df['is_valid'].sum()
invalid_smiles = total_smiles - valid_smiles
validity_pct = (valid_smiles / total_smiles) * 100

print(f"\n✅ Poprawne SMILES: {valid_smiles:,} ({validity_pct:.2f}%)")
print(f"❌ Niepoprawne SMILES: {invalid_smiles:,} ({(100-validity_pct):.2f}%)")

# Debug - pokaż przykłady
if valid_smiles > 0:
    print(f"\n✓ Przykłady POPRAWNYCH SMILES:")
    print(df[df['is_valid']]['SMILES'].head(3).tolist())
if invalid_smiles > 0:
    print(f"\n✗ Przykłady NIEPOPRAWNYCH SMILES:")
    print(df[~df['is_valid']]['SMILES'].head(3).tolist())

# ============================================================
# 2. ANALIZA BŁĘDÓW WALIDACJI
# ============================================================

if invalid_smiles > 0:
    print("\n" + "="*70)
    print("📊 ANALIZA BŁĘDÓW WALIDACJI")
    print("="*70)

    # Typy błędów
    error_counts = df[~df['is_valid']]['error_type'].value_counts()

    print("\n🔍 Typy błędów walidacji:")
    for error, count in error_counts.items():
        pct = (count / invalid_smiles) * 100
        print(f"  • {error:40s}  {count:5d} ({pct:5.2f}%)")

    # Przykłady błędnych SMILES
    print("\n⚠️  Przykłady niepoprawnych SMILES:")
    invalid_samples = df[~df['is_valid']].head(10)
    for idx, row in invalid_samples.iterrows():
        print(f"\n  [{idx}] {row['SMILES'][:60]}{'...' if len(row['SMILES']) > 60 else ''}")
        print(f"      Błąd: {row['error_type']}")

    # Wizualizacja typów błędów
    if len(error_counts) > 0:
        plt.figure(figsize=(12, 6))
        colors = plt.cm.Reds(np.linspace(0.4, 0.8, len(error_counts)))
        plt.barh(range(len(error_counts)), error_counts.values, color=colors, edgecolor='black')
        plt.yticks(range(len(error_counts)), error_counts.index)
        plt.xlabel('Liczba wystąpień', fontsize=12, fontweight='bold')
        plt.title('Typy błędów walidacji SMILES', fontsize=14, fontweight='bold')
        plt.grid(axis='x', alpha=0.3)

        for i, val in enumerate(error_counts.values):
            plt.text(val + max(error_counts.values)*0.01, i, f'{val}', 
                    va='center', fontweight='bold')

        plt.tight_layout()
        plt.savefig('qm9_smiles_errors.png', dpi=300, bbox_inches='tight')
        print("\n✅ Zapisano: qm9_smiles_errors.png")
        plt.show()

else:
    print("\n✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.")

# ============================================================
# 3. ANALIZA SKŁADU CHEMICZNEGO
# ============================================================

print("\n" + "="*70)
print("🧪 ANALIZA SKŁADU CHEMICZNEGO")
print("="*70)

# Tylko dla poprawnych SMILES
df_valid = df[df['is_valid']].copy()

def analyze_composition(smiles):
    """Analiza składu atomowego cząsteczki"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    atoms = [atom.GetSymbol() for atom in mol.GetAtoms()]
    return Counter(atoms)

print("\nAnalizuję skład atomowy cząsteczek...")
compositions = df_valid['SMILES'].apply(analyze_composition)

# Ekstrakcja wszystkich typów atomów
all_atoms = []
for comp in compositions.dropna():
    all_atoms.extend(comp.keys())

atom_types = Counter(all_atoms)

print(f"\n📊 Typy atomów występujące w datasecie ({len(atom_types)}):")
for atom, count in atom_types.most_common():
    pct = (count / len(df_valid)) * 100
    print(f"  • {atom:3s}  występuje w {count:6,} cząsteczkach ({pct:5.2f}%)")

# Rozkład liczby atomów per typ
atom_counts = {atom: [] for atom in atom_types.keys()}
for comp in compositions.dropna():
    for atom in atom_types.keys():
        atom_counts[atom].append(comp.get(atom, 0))

# Konwersja do DataFrame
df_atoms = pd.DataFrame(atom_counts)

print(f"\n📈 Statystyki liczby atomów:")
print(df_atoms.describe())

# ============================================================
# 4. WIZUALIZACJA SKŁADU ATOMOWEGO
# ============================================================

print("\n" + "="*70)
print("📊 WIZUALIZACJA SKŁADU ATOMOWEGO")
print("="*70)

# Wykresy dla głównych atomów (H, C, N, O, F)
main_atoms = ['H', 'C', 'N', 'O', 'F']
main_atoms = [atom for atom in main_atoms if atom in df_atoms.columns]

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()

for i, atom in enumerate(main_atoms):
    if atom in df_atoms.columns:
        data = df_atoms[atom]
        axes[i].hist(data, bins=50, color='steelblue', edgecolor='black', alpha=0.7)
        axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
        axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
        axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.1f}, Max: {data.max():.0f}', 
                         fontsize=12, fontweight='bold')
        axes[i].grid(axis='y', alpha=0.3)

# Wykres porównawczy - średnie liczby atomów
if len(main_atoms) < 6:
    ax_comparison = axes[5]
    means = [df_atoms[atom].mean() for atom in main_atoms]
    colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'][:len(main_atoms)]
    bars = ax_comparison.bar(main_atoms, means, color=colors, edgecolor='black', linewidth=1.5)
    ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
    ax_comparison.set_title('Średni skład atomowy cząsteczek QM9', fontsize=12, fontweight='bold')
    ax_comparison.grid(axis='y', alpha=0.3)

    for bar, val in zip(bars, means):
        ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.1, 
                          f'{val:.1f}', ha='center', fontweight='bold')

plt.tight_layout()
plt.savefig('qm9_atomic_composition.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_atomic_composition.png")
plt.show()

# ============================================================
# 5. ANALIZA WŁAŚCIWOŚCI MOLEKULARNYCH (RDKit)
# ============================================================

print("\n" + "="*70)
print("⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH")
print("="*70)

def calculate_rdkit_descriptors(smiles):
    """Oblicza podstawowe deskryptory RDKit"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    try:
        return {
            'MolWt': Descriptors.MolWt(mol),
            'NumAtoms': mol.GetNumAtoms(),
            'NumHeavyAtoms': mol.GetNumHeavyAtoms(),
            'NumHeteroatoms': rdMolDescriptors.CalcNumHeteroatoms(mol),
            'NumRotatableBonds': rdMolDescriptors.CalcNumRotatableBonds(mol),
            'NumHBD': rdMolDescriptors.CalcNumHBD(mol),
            'NumHBA': rdMolDescriptors.CalcNumHBA(mol),
            'NumRings': rdMolDescriptors.CalcNumRings(mol),
            'NumAromaticRings': rdMolDescriptors.CalcNumAromaticRings(mol),
            'TPSA': Descriptors.TPSA(mol),
            'LogP': Descriptors.MolLogP(mol),
            'FractionCSP3': rdMolDescriptors.CalcFractionCsp3(mol)
        }
    except:
        return None

print("Obliczam deskryptory RDKit dla poprawnych struktur...")
descriptors = df_valid['SMILES'].apply(calculate_rdkit_descriptors)
df_descriptors = descriptors.apply(pd.Series)

# Sprawdzenie czy są jakiekolwiek dane
if df_descriptors.empty or df_descriptors.shape[1] == 0:
    print("\n⚠️  BRAK POPRAWNYCH STRUKTUR DO ANALIZY!")
    print("   Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów.")
else:
    # Dodanie do głównego DataFrame
    df_valid = pd.concat([df_valid.reset_index(drop=True), 
                          df_descriptors.reset_index(drop=True)], axis=1)

    print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur")

    print("\n📊 Statystyki deskryptorów molekularnych:")
    print(df_descriptors.describe())

# ============================================================
# 6. CHARAKTERYSTYKA DATASETU QM9
# ============================================================

if not df_descriptors.empty and df_descriptors.shape[1] > 0:
    print("\n" + "="*70)
    print("📋 CHARAKTERYSTYKA DATASETU QM9")
    print("="*70)

    print(f"""
🔬 SKŁAD CHEMICZNY:
  • Dozwolone atomy: {', '.join(sorted(atom_types.keys()))}
  • Najczęstsze atomy: H, C, N, O, F (QM9 spec)
  • Średnia liczba atomów: {df_descriptors['NumAtoms'].mean():.1f}
  • Średnia liczba ciężkich atomów: {df_descriptors['NumHeavyAtoms'].mean():.1f}
  • Max ciężkich atomów: {df_descriptors['NumHeavyAtoms'].max():.0f} (limit QM9: 9)

📏 ROZMIAR MOLEKUŁ:
  • Masa molowa: {df_descriptors['MolWt'].mean():.1f} ± {df_descriptors['MolWt'].std():.1f} g/mol
  • Zakres: {df_descriptors['MolWt'].min():.1f} - {df_descriptors['MolWt'].max():.1f} g/mol

🔗 STRUKTURY:
  • Średnia liczba pierścieni: {df_descriptors['NumRings'].mean():.2f}
  • Cząsteczki aromatyczne: {(df_descriptors['NumAromaticRings'] > 0).sum():,} ({(df_descriptors['NumAromaticRings'] > 0).sum()/len(df_descriptors)*100:.1f}%)
  • Średnia rotowalne wiązania: {df_descriptors['NumRotatableBonds'].mean():.2f}

💧 WŁAŚCIWOŚCI FIZYKOCHEMICZNE:
  • Średnie LogP: {df_descriptors['LogP'].mean():.2f} ± {df_descriptors['LogP'].std():.2f}
  • Średnie TPSA: {df_descriptors['TPSA'].mean():.1f} Ų
  • Średnia frakcja sp3: {df_descriptors['FractionCSP3'].mean():.2f}

✅ JAKOŚĆ DANYCH:
  • Validacja SMILES: {validity_pct:.2f}%
  • Kompletność: {(df['is_valid'].sum() / len(df) * 100):.2f}%
  • Błędne struktury: {invalid_smiles:,}
""")

# ============================================================
# 7. WIZUALIZACJA PORÓWNAWCZA
# ============================================================

if not df_descriptors.empty and df_descriptors.shape[1] > 0:
    print("\n" + "="*70)
    print("📊 WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH")
    print("="*70)

    fig, axes = plt.subplots(2, 3, figsize=(18, 10))

    # 1. Masa molowa
    axes[0, 0].hist(df_descriptors['MolWt'], bins=50, color='#3498db', 
                    edgecolor='black', alpha=0.7)
    axes[0, 0].set_xlabel('Masa molowa (g/mol)', fontweight='bold')
    axes[0, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 0].set_title('Rozkład masy molowej', fontweight='bold', fontsize=12)
    axes[0, 0].grid(axis='y', alpha=0.3)

    # 2. Liczba ciężkich atomów
    axes[0, 1].hist(df_descriptors['NumHeavyAtoms'], bins=range(1, 11), 
                    color='#e74c3c', edgecolor='black', alpha=0.7)
    axes[0, 1].set_xlabel('Liczba ciężkich atomów', fontweight='bold')
    axes[0, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 1].set_title('Rozkład liczby ciężkich atomów (max=9)', fontweight='bold', fontsize=12)
    axes[0, 1].grid(axis='y', alpha=0.3)

    # 3. LogP
    axes[0, 2].hist(df_descriptors['LogP'], bins=50, color='#2ecc71', 
                    edgecolor='black', alpha=0.7)
    axes[0, 2].set_xlabel('LogP', fontweight='bold')
    axes[0, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 2].set_title('Rozkład lipofilowości (LogP)', fontweight='bold', fontsize=12)
    axes[0, 2].grid(axis='y', alpha=0.3)

    # 4. TPSA
    axes[1, 0].hist(df_descriptors['TPSA'], bins=50, color='#f39c12', 
                    edgecolor='black', alpha=0.7)
    axes[1, 0].set_xlabel('TPSA (Ų)', fontweight='bold')
    axes[1, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 0].set_title('Rozkład polarnej powierzchni (TPSA)', fontweight='bold', fontsize=12)
    axes[1, 0].grid(axis='y', alpha=0.3)

    # 5. Liczba pierścieni
    ring_counts = df_descriptors['NumRings'].value_counts().sort_index()
    axes[1, 1].bar(ring_counts.index, ring_counts.values, color='#9b59b6', 
                   edgecolor='black', alpha=0.7)
    axes[1, 1].set_xlabel('Liczba pierścieni', fontweight='bold')
    axes[1, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 1].set_title('Rozkład liczby pierścieni', fontweight='bold', fontsize=12)
    axes[1, 1].grid(axis='y', alpha=0.3)

    # 6. Frakcja sp3
    axes[1, 2].hist(df_descriptors['FractionCSP3'], bins=50, color='#1abc9c', 
                    edgecolor='black', alpha=0.7)
    axes[1, 2].set_xlabel('Frakcja węgli sp³', fontweight='bold')
    axes[1, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 2].set_title('Rozkład frakcji sp³ (saturacja)', fontweight='bold', fontsize=12)
    axes[1, 2].grid(axis='y', alpha=0.3)

    plt.suptitle('Właściwości molekularne datasetu QM9', fontsize=16, fontweight='bold', y=1.00)
    plt.tight_layout()
    plt.savefig('qm9_molecular_properties.png', dpi=300, bbox_inches='tight')
    print("✅ Zapisano: qm9_molecular_properties.png")
    plt.show()
else:
    print("\n⚠️  Brak danych do wizualizacji właściwości molekularnych")

print("\n" + "="*70)
print("✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!")
print("="*70)

======================================================================
WALIDACJA STRUKTUR SMILES - QM9 DATASET
======================================================================

📋 ETAP 1: PODSTAWOWA WALIDACJA
----------------------------------------------------------------------
Walidacja struktur SMILES...
Przykładowe SMILES do walidacji:
['[H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H]', '[H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H]', '[H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]1([H])OC([H])([H])[C@@]12C([H])([H])[H]', '[H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H]', '[H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2([H])[H]']

✅ Poprawne SMILES: 1,000 (100.00%)
❌ Niepoprawne SMILES: 0 (0.00%)

✓ Przykłady POPRAWNYCH SMILES:
['[H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H]', '[H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H]', '[H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]1([H])OC([H])([H])[C@@]12C([H])([H])[H]']

✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.

======================================================================
🧪 ANALIZA SKŁADU CHEMICZNEGO
======================================================================

Analizuję skład atomowy cząsteczek...

📊 Typy atomów występujące w datasecie (5):
  • C    występuje w  1,000 cząsteczkach (100.00%)
  • O    występuje w    843 cząsteczkach (84.30%)
  • N    występuje w    622 cząsteczkach (62.20%)
  • H    występuje w     16 cząsteczkach ( 1.60%)
  • F    występuje w     14 cząsteczkach ( 1.40%)

📈 Statystyki liczby atomów:
                 C            N            O            H            F
count  1000.000000  1000.000000  1000.000000  1000.000000  1000.000000
mean      6.251000     1.093000     1.392000     0.017000     0.020000
std       1.282046     1.138265     0.875844     0.136857     0.188774
min       2.000000     0.000000     0.000000     0.000000     0.000000
25%       5.000000     0.000000     1.000000     0.000000     0.000000
50%       6.000000     1.000000     1.000000     0.000000     0.000000
75%       7.000000     2.000000     2.000000     0.000000     0.000000
max       9.000000     6.000000     4.000000     2.000000     3.000000

======================================================================
📊 WIZUALIZACJA SKŁADU ATOMOWEGO
======================================================================
✅ Zapisano: qm9_atomic_composition.png

======================================================================
⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH
======================================================================
Obliczam deskryptory RDKit dla poprawnych struktur...

⚠️  BRAK POPRAWNYCH STRUKTUR DO ANALIZY!
   Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów.

⚠️  Brak danych do wizualizacji właściwości molekularnych

======================================================================
✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!
======================================================================

In [17]:

# ============================================================
# KOMÓRKA 7: WALIDACJA I ANALIZA STRUKTUR SMILES
# ============================================================

from rdkit import Chem, RDLogger
from rdkit.Chem import Descriptors, rdMolDescriptors
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter

# Wyciszenie ostrzeżeń RDKit
RDLogger.DisableLog('rdApp.*')

print("="*70)
print("WALIDACJA STRUKTUR SMILES - QM9 DATASET")
print("="*70)

# ============================================================
# 1. PODSTAWOWA WALIDACJA SMILES
# ============================================================

print("\n📋 ETAP 1: PODSTAWOWA WALIDACJA")
print("-" * 70)

# Usunięcie starych kolumn walidacji jeśli istnieją
if 'is_valid' in df.columns:
    df = df.drop(columns=['is_valid'])
if 'error_type' in df.columns:
    df = df.drop(columns=['error_type'])

def validate_smiles_detailed(smiles):
    """Szczegółowa walidacja SMILES z diagnozą błędów"""
    try:
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            return {'valid': False, 'error': 'Parse Error', 'mol': None}

        # Próba sanityzacji
        Chem.SanitizeMol(mol)
        return {'valid': True, 'error': None, 'mol': mol}

    except Exception as e:
        error_type = str(e).split(':')[0] if ':' in str(e) else str(e)
        return {'valid': False, 'error': error_type, 'mol': None}

# Walidacja wszystkich SMILES
print("Walidacja struktur SMILES...")
validation_results = df['SMILES'].apply(validate_smiles_detailed)

df['is_valid'] = validation_results.apply(lambda x: x['valid'])
df['error_type'] = validation_results.apply(lambda x: x['error'])

# Statystyki
total_smiles = len(df)
valid_smiles = df['is_valid'].sum()
invalid_smiles = total_smiles - valid_smiles
validity_pct = (valid_smiles / total_smiles) * 100

print(f"\n✅ Poprawne SMILES: {valid_smiles:,} ({validity_pct:.2f}%)")
print(f"❌ Niepoprawne SMILES: {invalid_smiles:,} ({(100-validity_pct):.2f}%)")

# ============================================================
# 2. SPRAWDZENIE CZY MAMY POPRAWNE DANE
# ============================================================

if valid_smiles == 0:
    print("\n" + "="*70)
    print("❌ BŁĄD: BRAK POPRAWNYCH STRUKTUR!")
    print("="*70)
    print("\nPrawdopodobne przyczyny:")
    print("  1. Kolumna 'is_valid' miała już wartości False z poprzedniego uruchomienia")
    print("  2. Problem z biblioteką RDKit")
    print("\nSpróbuj zrestartować kernel i uruchomić ponownie od początku.")
else:
    print(f"\n✓ Kontynuuję analizę {valid_smiles:,} poprawnych struktur...")

# Tylko dla poprawnych SMILES
df_valid = df[df['is_valid']].copy()

# ============================================================
# 3. ANALIZA SKŁADU CHEMICZNEGO
# ============================================================

if len(df_valid) > 0:
    print("\n" + "="*70)
    print("🧪 ANALIZA SKŁADU CHEMICZNEGO")
    print("="*70)

    def analyze_composition(smiles):
        """Analiza składu atomowego cząsteczki"""
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            return None
        atoms = [atom.GetSymbol() for atom in mol.GetAtoms()]
        return Counter(atoms)

    print("\nAnalizuję skład atomowy cząsteczek...")
    compositions = df_valid['SMILES'].apply(analyze_composition)

    # Ekstrakcja wszystkich typów atomów
    all_atoms = []
    for comp in compositions.dropna():
        all_atoms.extend(comp.keys())

    atom_types = Counter(all_atoms)

    print(f"\n📊 Typy atomów występujące w datasecie ({len(atom_types)}):")
    for atom, count in atom_types.most_common():
        pct = (count / len(df_valid)) * 100
        print(f"  • {atom:3s}  występuje w {count:6,} cząsteczkach ({pct:5.2f}%)")

    # Rozkład liczby atomów per typ
    atom_counts = {atom: [] for atom in atom_types.keys()}
    for comp in compositions.dropna():
        for atom in atom_types.keys():
            atom_counts[atom].append(comp.get(atom, 0))

    # Konwersja do DataFrame
    df_atoms = pd.DataFrame(atom_counts)

    print(f"\n📈 Statystyki liczby atomów:")
    print(df_atoms.describe())

    # ============================================================
    # 4. WIZUALIZACJA SKŁADU ATOMOWEGO
    # ============================================================

    print("\n" + "="*70)
    print("📊 WIZUALIZACJA SKŁADU ATOMOWEGO")
    print("="*70)

    # Wykresy dla głównych atomów (H, C, N, O, F)
    main_atoms = ['H', 'C', 'N', 'O', 'F']
    main_atoms = [atom for atom in main_atoms if atom in df_atoms.columns]

    fig, axes = plt.subplots(2, 3, figsize=(18, 10))
    axes = axes.flatten()

    for i, atom in enumerate(main_atoms):
        if atom in df_atoms.columns:
            data = df_atoms[atom]
            axes[i].hist(data, bins=50, color='steelblue', edgecolor='black', alpha=0.7)
            axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
            axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
            axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.1f}, Max: {data.max():.0f}', 
                             fontsize=12, fontweight='bold')
            axes[i].grid(axis='y', alpha=0.3)

    # Wykres porównawczy - średnie liczby atomów
    if len(main_atoms) < 6:
        ax_comparison = axes[5]
        means = [df_atoms[atom].mean() for atom in main_atoms]
        colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'][:len(main_atoms)]
        bars = ax_comparison.bar(main_atoms, means, color=colors, edgecolor='black', linewidth=1.5)
        ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
        ax_comparison.set_title('Średni skład atomowy cząsteczek QM9', fontsize=12, fontweight='bold')
        ax_comparison.grid(axis='y', alpha=0.3)

        for bar, val in zip(bars, means):
            ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.1, 
                              f'{val:.1f}', ha='center', fontweight='bold')

    plt.tight_layout()
    plt.savefig('qm9_atomic_composition.png', dpi=300, bbox_inches='tight')
    print("✅ Zapisano: qm9_atomic_composition.png")
    plt.show()

    # ============================================================
    # 5. ANALIZA WŁAŚCIWOŚCI MOLEKULARNYCH (RDKit)
    # ============================================================

    print("\n" + "="*70)
    print("⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH")
    print("="*70)

    def calculate_rdkit_descriptors(smiles):
        """Oblicza podstawowe deskryptory RDKit"""
        mol = Chem.MolFromSmiles(smiles)
        if mol is None:
            return None

        try:
            return {
                'MolWt': Descriptors.MolWt(mol),
                'NumAtoms': mol.GetNumAtoms(),
                'NumHeavyAtoms': mol.GetNumHeavyAtoms(),
                'NumHeteroatoms': rdMolDescriptors.CalcNumHeteroatoms(mol),
                'NumRotatableBonds': rdMolDescriptors.CalcNumRotatableBonds(mol),
                'NumHBD': rdMolDescriptors.CalcNumHBD(mol),
                'NumHBA': rdMolDescriptors.CalcNumHBA(mol),
                'NumRings': rdMolDescriptors.CalcNumRings(mol),
                'NumAromaticRings': rdMolDescriptors.CalcNumAromaticRings(mol),
                'TPSA': Descriptors.TPSA(mol),
                'LogP': Descriptors.MolLogP(mol),
                'FractionCSP3': rdMolDescriptors.CalcFractionCsp3(mol)
            }
        except:
            return None

    print("Obliczam deskryptory RDKit dla poprawnych struktur...")
    descriptors = df_valid['SMILES'].apply(calculate_rdkit_descriptors)
    df_descriptors = descriptors.apply(pd.Series)

    # Dodanie do głównego DataFrame
    df_valid = pd.concat([df_valid.reset_index(drop=True), 
                          df_descriptors.reset_index(drop=True)], axis=1)

    print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur")

    print("\n📊 Statystyki deskryptorów molekularnych:")
    print(df_descriptors.describe())

    # ============================================================
    # 6. WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH
    # ============================================================

    print("\n" + "="*70)
    print("📊 WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH")
    print("="*70)

    fig, axes = plt.subplots(2, 3, figsize=(18, 10))

    # 1. Masa molowa
    axes[0, 0].hist(df_descriptors['MolWt'], bins=50, color='#3498db', 
                    edgecolor='black', alpha=0.7)
    axes[0, 0].set_xlabel('Masa molowa (g/mol)', fontweight='bold')
    axes[0, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 0].set_title('Rozkład masy molowej', fontweight='bold', fontsize=12)
    axes[0, 0].grid(axis='y', alpha=0.3)

    # 2. Liczba ciężkich atomów
    axes[0, 1].hist(df_descriptors['NumHeavyAtoms'], bins=range(1, 11), 
                    color='#e74c3c', edgecolor='black', alpha=0.7)
    axes[0, 1].set_xlabel('Liczba ciężkich atomów', fontweight='bold')
    axes[0, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 1].set_title('Rozkład liczby ciężkich atomów (max=9)', fontweight='bold', fontsize=12)
    axes[0, 1].grid(axis='y', alpha=0.3)

    # 3. LogP
    axes[0, 2].hist(df_descriptors['LogP'], bins=50, color='#2ecc71', 
                    edgecolor='black', alpha=0.7)
    axes[0, 2].set_xlabel('LogP', fontweight='bold')
    axes[0, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[0, 2].set_title('Rozkład lipofilowości (LogP)', fontweight='bold', fontsize=12)
    axes[0, 2].grid(axis='y', alpha=0.3)

    # 4. TPSA
    axes[1, 0].hist(df_descriptors['TPSA'], bins=50, color='#f39c12', 
                    edgecolor='black', alpha=0.7)
    axes[1, 0].set_xlabel('TPSA (Ų)', fontweight='bold')
    axes[1, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 0].set_title('Rozkład polarnej powierzchni (TPSA)', fontweight='bold', fontsize=12)
    axes[1, 0].grid(axis='y', alpha=0.3)

    # 5. Liczba pierścieni
    ring_counts = df_descriptors['NumRings'].value_counts().sort_index()
    axes[1, 1].bar(ring_counts.index, ring_counts.values, color='#9b59b6', 
                   edgecolor='black', alpha=0.7)
    axes[1, 1].set_xlabel('Liczba pierścieni', fontweight='bold')
    axes[1, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 1].set_title('Rozkład liczby pierścieni', fontweight='bold', fontsize=12)
    axes[1, 1].grid(axis='y', alpha=0.3)

    # 6. Frakcja sp3
    axes[1, 2].hist(df_descriptors['FractionCSP3'], bins=50, color='#1abc9c', 
                    edgecolor='black', alpha=0.7)
    axes[1, 2].set_xlabel('Frakcja węgli sp³', fontweight='bold')
    axes[1, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
    axes[1, 2].set_title('Rozkład frakcji sp³ (saturacja)', fontweight='bold', fontsize=12)
    axes[1, 2].grid(axis='y', alpha=0.3)

    plt.suptitle('Właściwości molekularne datasetu QM9', fontsize=16, fontweight='bold', y=1.00)
    plt.tight_layout()
    plt.savefig('qm9_molecular_properties.png', dpi=300, bbox_inches='tight')
    print("✅ Zapisano: qm9_molecular_properties.png")
    plt.show()

print("\n" + "="*70)
print("✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!")
print("="*70)

======================================================================
WALIDACJA STRUKTUR SMILES - QM9 DATASET
======================================================================

📋 ETAP 1: PODSTAWOWA WALIDACJA
----------------------------------------------------------------------
Walidacja struktur SMILES...

✅ Poprawne SMILES: 1,000 (100.00%)
❌ Niepoprawne SMILES: 0 (0.00%)

✓ Kontynuuję analizę 1,000 poprawnych struktur...

======================================================================
🧪 ANALIZA SKŁADU CHEMICZNEGO
======================================================================

Analizuję skład atomowy cząsteczek...

📊 Typy atomów występujące w datasecie (5):
  • C    występuje w  1,000 cząsteczkach (100.00%)
  • O    występuje w    843 cząsteczkach (84.30%)
  • N    występuje w    622 cząsteczkach (62.20%)
  • H    występuje w     16 cząsteczkach ( 1.60%)
  • F    występuje w     14 cząsteczkach ( 1.40%)

📈 Statystyki liczby atomów:
                 C            N            O            H            F
count  1000.000000  1000.000000  1000.000000  1000.000000  1000.000000
mean      6.251000     1.093000     1.392000     0.017000     0.020000
std       1.282046     1.138265     0.875844     0.136857     0.188774
min       2.000000     0.000000     0.000000     0.000000     0.000000
25%       5.000000     0.000000     1.000000     0.000000     0.000000
50%       6.000000     1.000000     1.000000     0.000000     0.000000
75%       7.000000     2.000000     2.000000     0.000000     0.000000
max       9.000000     6.000000     4.000000     2.000000     3.000000

======================================================================
📊 WIZUALIZACJA SKŁADU ATOMOWEGO
======================================================================
✅ Zapisano: qm9_atomic_composition.png

======================================================================
⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH
======================================================================
Obliczam deskryptory RDKit dla poprawnych struktur...
✅ Obliczono 0 deskryptorów dla 1000 struktur

📊 Statystyki deskryptorów molekularnych:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[17], line 213
    210 print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur")
    212 print("\n📊 Statystyki deskryptorów molekularnych:")
--> 213 print(df_descriptors.describe())
    215 # ============================================================
    216 # 6. WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH
    217 # ============================================================
    219 print("\n" + "="*70)

File c:\Users\slast\miniconda3\envs\chem_eda\lib\site-packages\pandas\core\generic.py:12041, in NDFrame.describe(self, percentiles, include, exclude)
  11799 @final
  11800 def describe(
  11801     self,
   (...)
  11804     exclude=None,
  11805 ) -> Self:
  11806     """
  11807     Generate descriptive statistics.
  11808 
   (...)
  12039     max            NaN      3.0
  12040     """
> 12041     return describe_ndframe(
  12042         obj=self,
  12043         include=include,
  12044         exclude=exclude,
  12045         percentiles=percentiles,
  12046     ).__finalize__(self, method="describe")

File c:\Users\slast\miniconda3\envs\chem_eda\lib\site-packages\pandas\core\methods\describe.py:91, in describe_ndframe(obj, include, exclude, percentiles)
     87     describer = SeriesDescriber(
     88         obj=cast("Series", obj),
     89     )
     90 else:
---> 91     describer = DataFrameDescriber(
     92         obj=cast("DataFrame", obj),
     93         include=include,
     94         exclude=exclude,
     95     )
     97 result = describer.describe(percentiles=percentiles)
     98 return cast(NDFrameT, result)

File c:\Users\slast\miniconda3\envs\chem_eda\lib\site-packages\pandas\core\methods\describe.py:162, in DataFrameDescriber.__init__(self, obj, include, exclude)
    159 self.exclude = exclude
    161 if obj.ndim == 2 and obj.columns.size == 0:
--> 162     raise ValueError("Cannot describe a DataFrame without columns")
    164 super().__init__(obj)

ValueError: Cannot describe a DataFrame without columns

do naprawy potem¶

Rozkład train/valid/test - czy zbalansowany¶

In [19]:

# ============================================================
# KOMÓRKA 8: ANALIZA PODZIAŁU TRAIN/VALID/TEST
# ============================================================

print("="*70)
print("ANALIZA PODZIAŁU DATASETU - TRAIN/VALID/TEST")
print("="*70)

# ============================================================
# 1. PODSTAWOWE STATYSTYKI PODZIAŁU
# ============================================================

print("\n📊 ETAP 1: ROZKŁAD DANYCH WEDŁUG ZBIORÓW")
print("-" * 70)

# Liczebność poszczególnych zbiorów
split_counts = df['dataset_split'].value_counts().sort_index()
total_samples = len(df)

print("\n📈 Liczebność zbiorów:")
for split, count in split_counts.items():
    pct = (count / total_samples) * 100
    print(f"  • {split:6s}: {count:7,} próbek ({pct:5.2f}%)")

print(f"\n  RAZEM:    {total_samples:7,} próbek (100.00%)")

# ============================================================
# 2. WIZUALIZACJA PODZIAŁU
# ============================================================

print("\n" + "="*70)
print("📊 WIZUALIZACJA PODZIAŁU")
print("="*70)

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. Wykres kołowy
colors = ['#3498db', '#2ecc71', '#e74c3c']
wedges, texts, autotexts = axes[0].pie(split_counts.values, 
                                         labels=split_counts.index,
                                         autopct='%1.1f%%',
                                         colors=colors,
                                         startangle=90,
                                         textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[0].set_title('Rozkład procentowy zbiorów', fontsize=14, fontweight='bold')

# 2. Wykres słupkowy
bars = axes[1].bar(split_counts.index, split_counts.values, 
                    color=colors, edgecolor='black', linewidth=1.5, alpha=0.8)
axes[1].set_ylabel('Liczba próbek', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Zbiór', fontsize=12, fontweight='bold')
axes[1].set_title('Liczebność zbiorów', fontsize=14, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)

# Dodanie wartości na słupkach
for bar, val in zip(bars, split_counts.values):
    height = bar.get_height()
    axes[1].text(bar.get_x() + bar.get_width()/2., height,
                f'{val:,}', ha='center', va='bottom', fontweight='bold', fontsize=11)

# 3. Proporcje względem całości
proportions = (split_counts / total_samples * 100).values
axes[2].barh(split_counts.index, proportions, color=colors, 
             edgecolor='black', linewidth=1.5, alpha=0.8)
axes[2].set_xlabel('Procent całego datasetu (%)', fontsize=12, fontweight='bold')
axes[2].set_title('Proporcje zbiorów', fontsize=14, fontweight='bold')
axes[2].grid(axis='x', alpha=0.3)

# Dodanie wartości
for i, (split, val) in enumerate(zip(split_counts.index, proportions)):
    axes[2].text(val + 0.5, i, f'{val:.1f}%', 
                va='center', fontweight='bold', fontsize=11)

plt.tight_layout()
plt.savefig('qm9_dataset_split.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_dataset_split.png")
plt.show()

# ============================================================
# 3. ANALIZA BALANSU WŁAŚCIWOŚCI MIĘDZY ZBIORAMI
# ============================================================

print("\n" + "="*70)
print("🔬 ANALIZA BALANSU WŁAŚCIWOŚCI KWANTOWYCH")
print("="*70)

# Wybór kluczowych właściwości do analizy
key_properties = ['HOMO', 'LUMO', 'gap', 'dipole_moment', 'polarizability', 
                  'U0', 'H', 'G']

print("\n📊 Porównanie średnich wartości właściwości między zbiorami:")
print("-" * 70)

# Tabela porównawcza
comparison_data = []
for prop in key_properties:
    if prop in df.columns:
        train_mean = df[df['dataset_split'] == 'train'][prop].mean()
        valid_mean = df[df['dataset_split'] == 'valid'][prop].mean()
        test_mean = df[df['dataset_split'] == 'test'][prop].mean()

        # Obliczenie różnic procentowych
        max_diff = max(abs(train_mean - valid_mean), 
                       abs(train_mean - test_mean),
                       abs(valid_mean - test_mean))
        avg_mean = (train_mean + valid_mean + test_mean) / 3
        diff_pct = (max_diff / abs(avg_mean) * 100) if avg_mean != 0 else 0

        comparison_data.append({
            'Właściwość': prop,
            'Train': train_mean,
            'Valid': valid_mean,
            'Test': test_mean,
            'Max Diff %': diff_pct
        })

        print(f"\n{prop:25s}")
        print(f"  Train: {train_mean:10.4f}  |  Valid: {valid_mean:10.4f}  |  Test: {test_mean:10.4f}")
        print(f"  Max różnica: {diff_pct:.2f}%")

df_comparison = pd.DataFrame(comparison_data)

# Ocena balansu
print("\n" + "="*70)
print("📋 OCENA BALANSU ZBIORÓW")
print("="*70)

max_imbalance = df_comparison['Max Diff %'].max()

print(f"\n✓ Maksymalna różnica między zbiorami: {max_imbalance:.2f}%")

if max_imbalance < 1:
    balance_status = "🟢 DOSKONAŁY BALANS"
    balance_desc = "Zbiory są bardzo dobrze zbalansowane"
elif max_imbalance < 3:
    balance_status = "🟢 DOBRY BALANS"
    balance_desc = "Zbiory są dobrze zbalansowane"
elif max_imbalance < 5:
    balance_status = "🟡 AKCEPTOWALNY BALANS"
    balance_desc = "Niewielkie różnice między zbiorami"
else:
    balance_status = "🔴 SŁABY BALANS"
    balance_desc = "Znaczące różnice między zbiorami"

print(f"\n{balance_status}")
print(f"  {balance_desc}")

# ============================================================
# 4. WIZUALIZACJA ROZKŁADÓW WŁAŚCIWOŚCI
# ============================================================

print("\n" + "="*70)
print("📊 WIZUALIZACJA ROZKŁADÓW KLUCZOWYCH WŁAŚCIWOŚCI")
print("="*70)

# Wybór 6 najważniejszych właściwości
plot_properties = ['HOMO', 'LUMO', 'gap', 'dipole_moment', 'U0', 'H']

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()

for i, prop in enumerate(plot_properties):
    if prop in df.columns:
        # Rozkłady dla każdego zbioru
        for split, color in zip(['train', 'valid', 'test'], colors):
            data = df[df['dataset_split'] == split][prop].dropna()
            axes[i].hist(data, bins=50, alpha=0.5, label=split, 
                        color=color, edgecolor='black', linewidth=0.5)

        axes[i].set_xlabel(prop, fontsize=11, fontweight='bold')
        axes[i].set_ylabel('Liczba próbek', fontsize=11, fontweight='bold')
        axes[i].set_title(f'Rozkład {prop} w zbiorach', fontsize=12, fontweight='bold')
        axes[i].legend()
        axes[i].grid(axis='y', alpha=0.3)

plt.suptitle('Porównanie rozkładów właściwości między zbiorami', 
             fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_split_distributions.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_split_distributions.png")
plt.show()

# ============================================================
# 5. BOXPLOTY - PORÓWNANIE ZBIORÓW
# ============================================================

print("\n" + "="*70)
print("📦 BOXPLOTY - PORÓWNANIE MEDIANOWE")
print("="*70)

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()

for i, prop in enumerate(plot_properties):
    if prop in df.columns:
        # Przygotowanie danych
        data_to_plot = [df[df['dataset_split'] == split][prop].dropna() 
                        for split in ['train', 'valid', 'test']]

        bp = axes[i].boxplot(data_to_plot, 
                             labels=['Train', 'Valid', 'Test'],
                             patch_artist=True,
                             widths=0.6)

        # Kolorowanie boxów
        for patch, color in zip(bp['boxes'], colors):
            patch.set_facecolor(color)
            patch.set_alpha(0.7)

        axes[i].set_ylabel(prop, fontsize=11, fontweight='bold')
        axes[i].set_title(f'Boxplot {prop}', fontsize=12, fontweight='bold')
        axes[i].grid(axis='y', alpha=0.3)

plt.suptitle('Boxploty właściwości - porównanie zbiorów', 
             fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_split_boxplots.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_split_boxplots.png")
plt.show()

# ============================================================
# 6. PORÓWNANIE ODCHYLEŃ STANDARDOWYCH
# ============================================================

print("\n" + "="*70)
print("📊 PORÓWNANIE STATYSTYK OPISOWYCH")
print("="*70)

print("\nOdchylenia standardowe (im bardziej podobne, tym lepszy balans):")
print("-" * 70)

for prop in key_properties:
    if prop in df.columns:
        train_std = df[df['dataset_split'] == 'train'][prop].std()
        valid_std = df[df['dataset_split'] == 'valid'][prop].std()
        test_std = df[df['dataset_split'] == 'test'][prop].std()

        # Obliczenie różnicy
        max_std = max(train_std, valid_std, test_std)
        min_std = min(train_std, valid_std, test_std)
        diff_pct = ((max_std - min_std) / max_std * 100) if max_std > 0 else 0

        print(f"\n{prop:20s}")
        print(f"  Train: {train_std:10.4f}  |  Valid: {valid_std:10.4f}  |  Test: {test_std:10.4f}")
        print(f"  Różnica: {diff_pct:.2f}%")

# ============================================================
# 7. PODSUMOWANIE
# ============================================================

print("\n" + "="*70)
print("📋 PODSUMOWANIE ANALIZY PODZIAŁU")
print("="*70)

similar_count = (df_comparison['Max Diff %'] < 3).sum()
total_props = len(df_comparison)

print(f"""
📊 ROZKŁAD LICZBOWY:
  • Train: {split_counts['train']:,} ({split_counts['train']/total_samples*100:.1f}%)
  • Valid: {split_counts['valid']:,} ({split_counts['valid']/total_samples*100:.1f}%)
  • Test:  {split_counts['test']:,} ({split_counts['test']/total_samples*100:.1f}%)

⚖️  BALANS WŁAŚCIWOŚCI:
  • Status: {balance_status}
  • Maksymalna różnica średnich: {max_imbalance:.2f}%
  • Właściwości z dobrym balansem: {similar_count}/{total_props}

✅ WNIOSKI:
  • {'Podział jest dobrze zbalansowany pod względem liczebności' if max(split_counts.values)/min(split_counts.values) < 1.2 else 'Podział ma lekką nierówność liczebności'}
  • {'Rozkłady właściwości są spójne między zbiorami' if max_imbalance < 3 else 'Niektóre właściwości mogą różnić się między zbiorami'}
  • {'Dataset jest gotowy do trenowania modeli ML' if max_imbalance < 5 else 'Rozważ dodatkowe sprawdzenie podziału'}

💡 REKOMENDACJE:
  • {'Brak konieczności re-splittingu' if max_imbalance < 3 else 'Rozważ stratified split dla krytycznych właściwości'}
  • {'Możesz bezpiecznie trenować modele na tym podziale' if similar_count/total_props > 0.7 else 'Monitoruj metryki na wszystkich zbiorach podczas treningu'}
""")

print("="*70)
print("✅ ANALIZA PODZIAŁU ZAKOŃCZONA!")
print("="*70)

======================================================================
ANALIZA PODZIAŁU DATASETU - TRAIN/VALID/TEST
======================================================================

📊 ETAP 1: ROZKŁAD DANYCH WEDŁUG ZBIORÓW
----------------------------------------------------------------------

📈 Liczebność zbiorów:
  • test  :      91 próbek ( 9.10%)
  • train :     813 próbek (81.30%)
  • valid :      96 próbek ( 9.60%)

  RAZEM:      1,000 próbek (100.00%)

======================================================================
📊 WIZUALIZACJA PODZIAŁU
======================================================================
✅ Zapisano: qm9_dataset_split.png

======================================================================
🔬 ANALIZA BALANSU WŁAŚCIWOŚCI KWANTOWYCH
======================================================================

📊 Porównanie średnich wartości właściwości między zbiorami:
----------------------------------------------------------------------

HOMO                     
  Train:    -0.0224  |  Valid:    -0.0890  |  Test:    -0.0519
  Max różnica: 122.23%

LUMO                     
  Train:     0.0727  |  Valid:    -0.1290  |  Test:    -0.1003
  Max różnica: 386.37%

gap                      
  Train:     0.0839  |  Valid:    -0.0843  |  Test:    -0.0742
  Max różnica: 676.65%

dipole_moment            
  Train:    -0.0299  |  Valid:    -0.0693  |  Test:    -0.0120
  Max różnica: 154.55%

polarizability           
  Train:    -0.0297  |  Valid:    -0.0061  |  Test:    -0.0796
  Max różnica: 191.28%

U0                       
  Train:    -0.0125  |  Valid:    -0.0063  |  Test:    -0.0707
  Max różnica: 215.82%

H                        
  Train:     0.0289  |  Valid:     0.1364  |  Test:    -0.0878
  Max różnica: 867.44%

G                        
  Train:     0.0289  |  Valid:     0.1364  |  Test:    -0.0878
  Max różnica: 867.44%

======================================================================
📋 OCENA BALANSU ZBIORÓW
======================================================================

✓ Maksymalna różnica między zbiorami: 867.44%

🔴 SŁABY BALANS
  Znaczące różnice między zbiorami

======================================================================
📊 WIZUALIZACJA ROZKŁADÓW KLUCZOWYCH WŁAŚCIWOŚCI
======================================================================
✅ Zapisano: qm9_split_distributions.png

======================================================================
📦 BOXPLOTY - PORÓWNANIE MEDIANOWE
======================================================================

C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[i].boxplot(data_to_plot,
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[i].boxplot(data_to_plot,
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[i].boxplot(data_to_plot,
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[i].boxplot(data_to_plot,
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[i].boxplot(data_to_plot,
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[i].boxplot(data_to_plot,

✅ Zapisano: qm9_split_boxplots.png

======================================================================
📊 PORÓWNANIE STATYSTYK OPISOWYCH
======================================================================

Odchylenia standardowe (im bardziej podobne, tym lepszy balans):
----------------------------------------------------------------------

HOMO                
  Train:     1.0495  |  Valid:     0.9226  |  Test:     0.9963
  Różnica: 12.09%

LUMO                
  Train:     0.9760  |  Valid:     1.0129  |  Test:     1.1192
  Różnica: 12.80%

gap                 
  Train:     1.0070  |  Valid:     0.9814  |  Test:     1.0664
  Różnica: 7.97%

dipole_moment       
  Train:     1.0284  |  Valid:     0.9200  |  Test:     0.9418
  Różnica: 10.55%

polarizability      
  Train:     1.0579  |  Valid:     1.1738  |  Test:     1.0389
  Różnica: 11.50%

U0                  
  Train:     1.0377  |  Valid:     1.1266  |  Test:     1.0621
  Różnica: 7.89%

H                   
  Train:     1.0520  |  Valid:     1.1977  |  Test:     1.0812
  Różnica: 12.17%

G                   
  Train:     1.0520  |  Valid:     1.1977  |  Test:     1.0812
  Różnica: 12.17%

======================================================================
📋 PODSUMOWANIE ANALIZY PODZIAŁU
======================================================================

📊 ROZKŁAD LICZBOWY:
  • Train: 813 (81.3%)
  • Valid: 96 (9.6%)
  • Test:  91 (9.1%)

⚖️  BALANS WŁAŚCIWOŚCI:
  • Status: 🔴 SŁABY BALANS
  • Maksymalna różnica średnich: 867.44%
  • Właściwości z dobrym balansem: 0/8

✅ WNIOSKI:
  • Podział ma lekką nierówność liczebności
  • Niektóre właściwości mogą różnić się między zbiorami
  • Rozważ dodatkowe sprawdzenie podziału

💡 REKOMENDACJE:
  • Rozważ stratified split dla krytycznych właściwości
  • Monitoruj metryki na wszystkich zbiorach podczas treningu

======================================================================
✅ ANALIZA PODZIAŁU ZAKOŃCZONA!
======================================================================

3. ANALIZA SKŁADU CHEMICZNEGO 🧪¶

• Rozkład atomów (H, C, N, O, F) - ile każdego typu
• Masa molowa - rozkład, median, percentyle
• Liczba atomów ciężkich - histogram, statystyki
• Liczba atomów wodoru - analiza
• Wzory sumaryczne - najczęstsze składy

In [20]:

# ============================================================
# KOMÓRKA 9: ANALIZA SKŁADU CHEMICZNEGO
# ============================================================

print("="*70)
print("ANALIZA SKŁADU CHEMICZNEGO - QM9 DATASET")
print("="*70)

# ============================================================
# 1. ROZKŁAD ATOMÓW (H, C, N, O, F)
# ============================================================

print("\n🧪 ETAP 1: ROZKŁAD ATOMÓW W DATASECIE")
print("-" * 70)

# Funkcja do analizy składu
def get_atom_composition(smiles):
    """Zwraca słownik z liczbą atomów każdego typu"""
    from rdkit import Chem
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    atom_dict = {'H': 0, 'C': 0, 'N': 0, 'O': 0, 'F': 0}
    for atom in mol.GetAtoms():
        symbol = atom.GetSymbol()
        if symbol in atom_dict:
            atom_dict[symbol] += 1

    # Dodaj wodory (mogą być jawne lub niejawne)
    atom_dict['H'] = sum(1 for atom in mol.GetAtoms() if atom.GetSymbol() == 'H')
    if atom_dict['H'] == 0:  # Jeśli wodory są niejawne
        atom_dict['H'] = sum(atom.GetTotalNumHs() for atom in mol.GetAtoms())

    return atom_dict

# Obliczenie składu dla wszystkich cząsteczek
print("Analizuję skład atomowy cząsteczek...")
compositions = df['SMILES'].apply(get_atom_composition)
df_composition = compositions.apply(pd.Series)

# Dodanie do głównego DataFrame
df['H_count'] = df_composition['H']
df['C_count'] = df_composition['C']
df['N_count'] = df_composition['N']
df['O_count'] = df_composition['O']
df['F_count'] = df_composition['F']

print("\n📊 Statystyki rozkładu atomów:")
print("-" * 70)

atom_stats = df_composition.describe()
print(atom_stats)

# Całkowita liczba każdego atomu w datasecie
print("\n📈 Całkowita liczba atomów w datasecie:")
for atom in ['H', 'C', 'N', 'O', 'F']:
    total = df[f'{atom}_count'].sum()
    avg_per_mol = df[f'{atom}_count'].mean()
    print(f"  • {atom}: {total:,} atomów | Średnio {avg_per_mol:.2f} na cząsteczkę")

# ============================================================
# 2. WIZUALIZACJA ROZKŁADU ATOMÓW
# ============================================================

print("\n" + "="*70)
print("📊 WIZUALIZACJA ROZKŁADU ATOMÓW")
print("="*70)

fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()

atom_colors = {
    'H': '#E8E8E8',  # Jasny szary (wodór)
    'C': '#909090',  # Ciemny szary (węgiel)
    'N': '#3050F8',  # Niebieski (azot)
    'O': '#FF0D0D',  # Czerwony (tlen)
    'F': '#90E050'   # Zielony (fluor)
}

for i, atom in enumerate(['H', 'C', 'N', 'O', 'F']):
    col = f'{atom}_count'
    data = df[col].dropna()

    axes[i].hist(data, bins=50, color=atom_colors[atom], 
                 edgecolor='black', alpha=0.7, linewidth=0.5)
    axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
    axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
    axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.2f} | Max: {data.max():.0f}', 
                      fontsize=12, fontweight='bold')
    axes[i].grid(axis='y', alpha=0.3)

    # Dodanie linii mediany
    median_val = data.median()
    axes[i].axvline(median_val, color='red', linestyle='--', linewidth=2, 
                    label=f'Mediana: {median_val:.1f}')
    axes[i].legend()

# Wykres porównawczy - średnie liczby atomów
ax_comparison = axes[5]
atoms = ['H', 'C', 'N', 'O', 'F']
means = [df[f'{atom}_count'].mean() for atom in atoms]
colors_list = [atom_colors[atom] for atom in atoms]

bars = ax_comparison.bar(atoms, means, color=colors_list, 
                          edgecolor='black', linewidth=1.5, alpha=0.8)
ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
ax_comparison.set_xlabel('Typ atomu', fontsize=11, fontweight='bold')
ax_comparison.set_title('Średni skład atomowy cząsteczek QM9', 
                        fontsize=12, fontweight='bold')
ax_comparison.grid(axis='y', alpha=0.3)

for bar, val in zip(bars, means):
    ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.2, 
                       f'{val:.2f}', ha='center', fontweight='bold', fontsize=11)

plt.suptitle('Rozkład atomów w datasecie QM9', fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_atom_distribution.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_atom_distribution.png")
plt.show()

# ============================================================
# 3. MASA MOLOWA - SZCZEGÓŁOWA ANALIZA
# ============================================================

print("\n" + "="*70)
print("⚖️  ANALIZA MASY MOLOWEJ")
print("="*70)

# Obliczenie masy molowej
from rdkit.Chem import Descriptors

def calculate_molecular_weight(smiles):
    from rdkit import Chem
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    return Descriptors.MolWt(mol)

df['MolWt'] = df['SMILES'].apply(calculate_molecular_weight)

# Statystyki
print("\n📊 Statystyki masy molowej:")
print(f"  • Średnia:    {df['MolWt'].mean():.2f} g/mol")
print(f"  • Mediana:    {df['MolWt'].median():.2f} g/mol")
print(f"  • Odch. std:  {df['MolWt'].std():.2f} g/mol")
print(f"  • Min:        {df['MolWt'].min():.2f} g/mol")
print(f"  • Max:        {df['MolWt'].max():.2f} g/mol")

# Percentyle
print("\n📈 Percentyle masy molowej:")
percentiles = [10, 25, 50, 75, 90, 95, 99]
for p in percentiles:
    val = df['MolWt'].quantile(p/100)
    print(f"  • {p:2d}%:  {val:.2f} g/mol")

# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. Histogram
axes[0].hist(df['MolWt'], bins=60, color='steelblue', 
             edgecolor='black', alpha=0.7, linewidth=0.5)
axes[0].axvline(df['MolWt'].mean(), color='red', linestyle='--', 
                linewidth=2, label=f"Średnia: {df['MolWt'].mean():.1f}")
axes[0].axvline(df['MolWt'].median(), color='green', linestyle='--', 
                linewidth=2, label=f"Mediana: {df['MolWt'].median():.1f}")
axes[0].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład masy molowej', fontsize=13, fontweight='bold')
axes[0].legend()
axes[0].grid(axis='y', alpha=0.3)

# 2. Boxplot
bp = axes[1].boxplot([df['MolWt'].dropna()], labels=['Masa molowa'], 
                      patch_artist=True, widths=0.5)
bp['boxes'][0].set_facecolor('lightblue')
bp['boxes'][0].set_alpha(0.7)
axes[1].set_ylabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[1].set_title('Boxplot masy molowej', fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)

# 3. Dystrybuanta (CDF)
sorted_mw = np.sort(df['MolWt'].dropna())
cdf = np.arange(1, len(sorted_mw) + 1) / len(sorted_mw)
axes[2].plot(sorted_mw, cdf, linewidth=2, color='darkblue')
axes[2].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[2].set_ylabel('Dystrybuanta (CDF)', fontsize=12, fontweight='bold')
axes[2].set_title('Dystrybuanta masy molowej', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)

# Dodanie linii percentyli
for p in [25, 50, 75]:
    val = df['MolWt'].quantile(p/100)
    axes[2].axvline(val, color='red', linestyle=':', alpha=0.5)
    axes[2].text(val, 0.05, f'{p}%', rotation=90, fontsize=9)

plt.tight_layout()
plt.savefig('qm9_molecular_weight.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_molecular_weight.png")
plt.show()

# ============================================================
# 4. LICZBA ATOMÓW CIĘŻKICH
# ============================================================

print("\n" + "="*70)
print("🔬 ANALIZA ATOMÓW CIĘŻKICH (non-H)")
print("="*70)

# Obliczenie liczby atomów ciężkich
df['HeavyAtoms'] = df['C_count'] + df['N_count'] + df['O_count'] + df['F_count']

print("\n📊 Statystyki atomów ciężkich:")
print(f"  • Średnia:   {df['HeavyAtoms'].mean():.2f} atomów")
print(f"  • Mediana:   {df['HeavyAtoms'].median():.0f} atomów")
print(f"  • Min:       {df['HeavyAtoms'].min():.0f} atomów")
print(f"  • Max:       {df['HeavyAtoms'].max():.0f} atomów")
print(f"  • Limit QM9: 9 atomów ciężkich")

# Rozkład
print("\n📈 Rozkład liczby atomów ciężkich:")
heavy_dist = df['HeavyAtoms'].value_counts().sort_index()
for count, freq in heavy_dist.items():
    pct = (freq / len(df)) * 100
    bar = '█' * int(pct / 2)
    print(f"  {count:2.0f} atomów: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")

# Wizualizacja
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Histogram
axes[0].bar(heavy_dist.index, heavy_dist.values, color='coral', 
            edgecolor='black', alpha=0.8, linewidth=1)
axes[0].set_xlabel('Liczba atomów ciężkich', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład liczby atomów ciężkich\n(QM9: max 9 atomów)', 
                  fontsize=13, fontweight='bold')
axes[0].grid(axis='y', alpha=0.3)
axes[0].axvline(9, color='red', linestyle='--', linewidth=2, label='Limit QM9: 9')
axes[0].legend()

# Pie chart dla grup
bins = [0, 3, 6, 9, 12]
labels = ['1-3 atomy', '4-6 atomów', '7-9 atomów', '10+ atomów']
df['HeavyAtoms_bin'] = pd.cut(df['HeavyAtoms'], bins=bins, labels=labels[:len(bins)-1])
heavy_grouped = df['HeavyAtoms_bin'].value_counts()

axes[1].pie(heavy_grouped.values, labels=heavy_grouped.index, autopct='%1.1f%%',
            colors=['#ff9999', '#66b3ff', '#99ff99'], startangle=90,
            textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[1].set_title('Grupy wielkości molekuł', fontsize=13, fontweight='bold')

plt.tight_layout()
plt.savefig('qm9_heavy_atoms.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_heavy_atoms.png")
plt.show()

# ============================================================
# 5. WZORY SUMARYCZNE - NAJCZĘSTSZE SKŁADY
# ============================================================

print("\n" + "="*70)
print("🧬 WZORY SUMARYCZNE - NAJCZĘSTSZE SKŁADY")
print("="*70)

# Tworzenie wzorów sumarycznych
def create_formula(row):
    """Tworzy wzór sumaryczny w formacie CxHyNzOaFb"""
    parts = []
    for atom in ['C', 'H', 'N', 'O', 'F']:
        count = int(row[f'{atom}_count'])
        if count > 0:
            if count == 1:
                parts.append(atom)
            else:
                parts.append(f"{atom}{count}")
    return ''.join(parts) if parts else "Empty"

df['MolecularFormula'] = df.apply(create_formula, axis=1)

# Top 20 najczęstszych wzorów
print("\n🏆 TOP 20 najczęstszych wzorów sumarycznych:")
print("-" * 70)

top_formulas = df['MolecularFormula'].value_counts().head(20)

for i, (formula, count) in enumerate(top_formulas.items(), 1):
    pct = (count / len(df)) * 100
    print(f"  {i:2d}. {formula:15s}  {count:4,} cząsteczek ({pct:5.2f}%)")

# Wizualizacja top 15
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# Barplot
top15 = top_formulas.head(15)
axes[0].barh(range(len(top15)), top15.values, color='teal', 
             edgecolor='black', alpha=0.7)
axes[0].set_yticks(range(len(top15)))
axes[0].set_yticklabels(top15.index, fontsize=9)
axes[0].set_xlabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Top 15 najczęstszych wzorów sumarycznych', 
                  fontsize=13, fontweight='bold')
axes[0].grid(axis='x', alpha=0.3)
axes[0].invert_yaxis()

# Dodanie wartości
for i, val in enumerate(top15.values):
    axes[0].text(val + 5, i, f'{val}', va='center', fontsize=9, fontweight='bold')

# Rozkład unikalności
unique_formulas = df['MolecularFormula'].nunique()
total_molecules = len(df)

categories = ['Unikalne wzory', 'Powtarzające się']
values = [unique_formulas, total_molecules - unique_formulas]

axes[1].pie(values, labels=categories, autopct='%1.1f%%',
            colors=['#ff6b6b', '#4ecdc4'], startangle=90,
            textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[1].set_title(f'Unikalność wzorów sumarycznych\n({unique_formulas:,} unikalnych / {total_molecules:,} total)', 
                  fontsize=13, fontweight='bold')

plt.tight_layout()
plt.savefig('qm9_molecular_formulas.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_molecular_formulas.png")
plt.show()

# ============================================================
# 6. PODSUMOWANIE SKŁADU CHEMICZNEGO
# ============================================================

print("\n" + "="*70)
print("📋 PODSUMOWANIE SKŁADU CHEMICZNEGO")
print("="*70)

print(f"""
🧪 SKŁAD ATOMOWY:
  • Średnia H:  {df['H_count'].mean():.2f} (wodór)
  • Średnia C:  {df['C_count'].mean():.2f} (węgiel)
  • Średnia N:  {df['N_count'].mean():.2f} (azot)
  • Średnia O:  {df['O_count'].mean():.2f} (tlen)
  • Średnia F:  {df['F_count'].mean():.2f} (fluor)

⚖️  MASA MOLOWA:
  • Zakres:     {df['MolWt'].min():.1f} - {df['MolWt'].max():.1f} g/mol
  • Średnia:    {df['MolWt'].mean():.1f} g/mol
  • Mediana:    {df['MolWt'].median():.1f} g/mol

🔬 ATOMY CIĘŻKIE:
  • Zakres:     {df['HeavyAtoms'].min():.0f} - {df['HeavyAtoms'].max():.0f} atomów
  • Średnia:    {df['HeavyAtoms'].mean():.2f} atomów
  • Limit QM9:  9 atomów ciężkich (spec datasetu)

🧬 WZORY SUMARYCZNE:
  • Unikalne wzory:           {unique_formulas:,}
  • Całkowita liczba molekuł: {total_molecules:,}
  • Najczęstszy wzór:         {top_formulas.index[0]} ({top_formulas.values[0]} molekuł)

✅ CHARAKTERYSTYKA DATASETU:
  • QM9 zawiera małe cząsteczki organiczne (max 9 ciężkich atomów)
  • Dominują związki z H, C, N, O, F (jak w specyfikacji)
  • Duża różnorodność strukturalna ({unique_formulas:,} unikalnych wzorów)
  • Masa molowa typowo poniżej 150 g/mol (small molecules)
""")

print("="*70)
print("✅ ANALIZA SKŁADU CHEMICZNEGO ZAKOŃCZONA!")
print("="*70)

======================================================================
ANALIZA SKŁADU CHEMICZNEGO - QM9 DATASET
======================================================================

🧪 ETAP 1: ROZKŁAD ATOMÓW W DATASECIE
----------------------------------------------------------------------
Analizuję skład atomowy cząsteczek...

📊 Statystyki rozkładu atomów:
----------------------------------------------------------------------
                 H            C            N            O            F
count  1000.000000  1000.000000  1000.000000  1000.000000  1000.000000
mean      9.322000     6.251000     1.093000     1.392000     0.020000
std       3.103215     1.282046     1.138265     0.875844     0.188774
min       1.000000     2.000000     0.000000     0.000000     0.000000
25%       7.000000     5.000000     0.000000     1.000000     0.000000
50%       9.000000     6.000000     1.000000     1.000000     0.000000
75%      12.000000     7.000000     2.000000     2.000000     0.000000
max      18.000000     9.000000     6.000000     4.000000     3.000000

📈 Całkowita liczba atomów w datasecie:
  • H: 9,322 atomów | Średnio 9.32 na cząsteczkę
  • C: 6,251 atomów | Średnio 6.25 na cząsteczkę
  • N: 1,093 atomów | Średnio 1.09 na cząsteczkę
  • O: 1,392 atomów | Średnio 1.39 na cząsteczkę
  • F: 20 atomów | Średnio 0.02 na cząsteczkę

======================================================================
📊 WIZUALIZACJA ROZKŁADU ATOMÓW
======================================================================
✅ Zapisano: qm9_atom_distribution.png

======================================================================
⚖️  ANALIZA MASY MOLOWEJ
======================================================================

📊 Statystyki masy molowej:
  • Średnia:    122.53 g/mol
  • Mediana:    125.13 g/mol
  • Odch. std:  8.67 g/mol
  • Min:        59.07 g/mol
  • Max:        142.12 g/mol

📈 Percentyle masy molowej:
  • 10%:  112.11 g/mol
  • 25%:  121.10 g/mol
  • 50%:  125.13 g/mol
  • 75%:  128.13 g/mol
  • 90%:  130.08 g/mol
  • 95%:  131.13 g/mol
  • 99%:  132.20 g/mol

C:\Users\slast\AppData\Local\Temp\ipykernel_6436\3117465052.py:175: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[1].boxplot([df['MolWt'].dropna()], labels=['Masa molowa'],

✅ Zapisano: qm9_molecular_weight.png

======================================================================
🔬 ANALIZA ATOMÓW CIĘŻKICH (non-H)
======================================================================

📊 Statystyki atomów ciężkich:
  • Średnia:   8.76 atomów
  • Mediana:   9 atomów
  • Min:       4 atomów
  • Max:       9 atomów
  • Limit QM9: 9 atomów ciężkich

📈 Rozkład liczby atomów ciężkich:
   4 atomów:      1 cząsteczek ( 0.10%) 
   5 atomów:      2 cząsteczek ( 0.20%) 
   6 atomów:      8 cząsteczek ( 0.80%) 
   7 atomów:     30 cząsteczek ( 3.00%) █
   8 atomów:    147 cząsteczek (14.70%) ███████
   9 atomów:    812 cząsteczek (81.20%) ████████████████████████████████████████

✅ Zapisano: qm9_heavy_atoms.png

======================================================================
🧬 WZORY SUMARYCZNE - NAJCZĘSTSZE SKŁADY
======================================================================

🏆 TOP 20 najczęstszych wzorów sumarycznych:
----------------------------------------------------------------------
   1. C7H10O2            52 cząsteczek ( 5.20%)
   2. C6H9NO2            39 cząsteczek ( 3.90%)
   3. C7H12O2            36 cząsteczek ( 3.60%)
   4. C8H12O             35 cząsteczek ( 3.50%)
   5. C7H9NO             25 cząsteczek ( 2.50%)
   6. C6H8O3             23 cząsteczek ( 2.30%)
   7. C8H10O             21 cząsteczek ( 2.10%)
   8. C6H7NO2            19 cząsteczek ( 1.90%)
   9. C7H8O2             19 cząsteczek ( 1.90%)
  10. C7H14NO            18 cząsteczek ( 1.80%)
  11. C6H10N2O           17 cząsteczek ( 1.70%)
  12. C6H10O3            16 cząsteczek ( 1.60%)
  13. C7H12NO            15 cząsteczek ( 1.50%)
  14. C8H14O             14 cząsteczek ( 1.40%)
  15. C6H8N2O            14 cząsteczek ( 1.40%)
  16. C6H11NO2           13 cząsteczek ( 1.30%)
  17. C6H12O3            13 cząsteczek ( 1.30%)
  18. C8H16O             12 cząsteczek ( 1.20%)
  19. C7H7NO             12 cząsteczek ( 1.20%)
  20. C9H14              11 cząsteczek ( 1.10%)

✅ Zapisano: qm9_molecular_formulas.png

======================================================================
📋 PODSUMOWANIE SKŁADU CHEMICZNEGO
======================================================================

🧪 SKŁAD ATOMOWY:
  • Średnia H:  9.32 (wodór)
  • Średnia C:  6.25 (węgiel)
  • Średnia N:  1.09 (azot)
  • Średnia O:  1.39 (tlen)
  • Średnia F:  0.02 (fluor)

⚖️  MASA MOLOWA:
  • Zakres:     59.1 - 142.1 g/mol
  • Średnia:    122.5 g/mol
  • Mediana:    125.1 g/mol

🔬 ATOMY CIĘŻKIE:
  • Zakres:     4 - 9 atomów
  • Średnia:    8.76 atomów
  • Limit QM9:  9 atomów ciężkich (spec datasetu)

🧬 WZORY SUMARYCZNE:
  • Unikalne wzory:           243
  • Całkowita liczba molekuł: 1,000
  • Najczęstszy wzór:         C7H10O2 (52 molekuł)

✅ CHARAKTERYSTYKA DATASETU:
  • QM9 zawiera małe cząsteczki organiczne (max 9 ciężkich atomów)
  • Dominują związki z H, C, N, O, F (jak w specyfikacji)
  • Duża różnorodność strukturalna (243 unikalnych wzorów)
  • Masa molowa typowo poniżej 150 g/mol (small molecules)

======================================================================
✅ ANALIZA SKŁADU CHEMICZNEGO ZAKOŃCZONA!
======================================================================

4. ANALIZA STRUKTURALNA 🔬¶

• Liczba wiązań - single, double, triple, aromatic
• Liczba pierścieni - histogram, rozkład 0-3 pierścienie
• Pierścienie aromatyczne vs alifatyczne
• Rotowalne wiązania - flexibility molekuł
• Frakcja sp3 - saturacja vs aromatyczność
• Chiralność - czy są centra chiralne

In [25]:

# ============================================================
# KOMÓRKA 10: ANALIZA STRUKTURALNA
# ============================================================

print("="*70)
print("ANALIZA STRUKTURALNA - QM9 DATASET")
print("="*70)

# ============================================================
# 1. ANALIZA WIĄZAŃ (SINGLE, DOUBLE, TRIPLE, AROMATIC)
# ============================================================

print("\n🔗 ETAP 1: ANALIZA TYPÓW WIĄZAŃ")
print("-" * 70)

from rdkit import Chem
from rdkit.Chem import rdMolDescriptors, Descriptors

def analyze_bonds(smiles):
    """Analiza typów wiązań w cząsteczce"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    bond_types = {
        'single': 0,
        'double': 0,
        'triple': 0,
        'aromatic': 0
    }

    for bond in mol.GetBonds():
        bond_type = bond.GetBondType()
        if bond_type == Chem.BondType.SINGLE:
            if not bond.GetIsAromatic():
                bond_types['single'] += 1
        elif bond_type == Chem.BondType.DOUBLE:
            if not bond.GetIsAromatic():
                bond_types['double'] += 1
        elif bond_type == Chem.BondType.TRIPLE:
            bond_types['triple'] += 1

        if bond.GetIsAromatic():
            bond_types['aromatic'] += 1

    return bond_types

print("Analizuję typy wiązań w cząsteczkach...")
bonds_data = df['SMILES'].apply(analyze_bonds)
df_bonds = bonds_data.apply(pd.Series)

# Dodanie do głównego DataFrame
df['single_bonds'] = df_bonds['single']
df['double_bonds'] = df_bonds['double']
df['triple_bonds'] = df_bonds['triple']
df['aromatic_bonds'] = df_bonds['aromatic']

print("\n📊 Statystyki typów wiązań:")
print(df_bonds.describe())

print("\n📈 Średnia liczba wiązań każdego typu:")
for bond_type in ['single', 'double', 'triple', 'aromatic']:
    avg = df[f'{bond_type}_bonds'].mean()
    total = df[f'{bond_type}_bonds'].sum()
    molecules_with = (df[f'{bond_type}_bonds'] > 0).sum()
    pct = (molecules_with / len(df)) * 100
    print(f"  • {bond_type:10s}: Średnio {avg:.2f} | Molekuły z tym wiązaniem: {molecules_with:,} ({pct:.1f}%)")

# Wizualizacja
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

bond_colors = {
    'single': '#3498db',
    'double': '#e74c3c',
    'triple': '#f39c12',
    'aromatic': '#9b59b6'
}

for idx, (bond_type, ax) in enumerate(zip(['single', 'double', 'triple', 'aromatic'], axes.flatten())):
    data = df[f'{bond_type}_bonds'].dropna()
    ax.hist(data, bins=40, color=bond_colors[bond_type], 
            edgecolor='black', alpha=0.7, linewidth=0.5)
    ax.set_xlabel(f'Liczba wiązań {bond_type}', fontsize=11, fontweight='bold')
    ax.set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
    ax.set_title(f'Rozkład wiązań {bond_type}\nŚrednia: {data.mean():.2f}', 
                 fontsize=12, fontweight='bold')
    ax.grid(axis='y', alpha=0.3)
    ax.axvline(data.mean(), color='red', linestyle='--', linewidth=2, label='Średnia')
    ax.legend()

plt.suptitle('Rozkład typów wiązań w QM9', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_bond_types.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_bond_types.png")
plt.show()

# ============================================================
# 2. LICZBA PIERŚCIENI
# ============================================================

print("\n" + "="*70)
print("⭕ ANALIZA PIERŚCIENI")
print("="*70)

def analyze_rings(smiles):
    """Szczegółowa analiza pierścieni"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    return {
        'total_rings': rdMolDescriptors.CalcNumRings(mol),
        'aromatic_rings': rdMolDescriptors.CalcNumAromaticRings(mol),
        'aliphatic_rings': rdMolDescriptors.CalcNumAliphaticRings(mol),
        'saturated_rings': rdMolDescriptors.CalcNumSaturatedRings(mol),
        'hetero_rings': rdMolDescriptors.CalcNumHeterocycles(mol),
        'aromatic_hetero_rings': rdMolDescriptors.CalcNumAromaticHeterocycles(mol)
    }

print("Analizuję struktury pierścieniowe...")
rings_data = df['SMILES'].apply(analyze_rings)
df_rings = rings_data.apply(pd.Series)

# Dodanie do głównego DataFrame
for col in df_rings.columns:
    df[col] = df_rings[col]

print("\n📊 Statystyki pierścieni:")
print(df_rings.describe())

# Rozkład liczby pierścieni
print("\n📈 Rozkład liczby pierścieni:")
ring_dist = df['total_rings'].value_counts().sort_index()
for count, freq in ring_dist.items():
    pct = (freq / len(df)) * 100
    bar = '█' * int(pct / 3)
    print(f"  {count:1.0f} pierścieni: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")

# Molekuły z/bez pierścieni
molecules_with_rings = (df['total_rings'] > 0).sum()
molecules_without_rings = (df['total_rings'] == 0).sum()
print(f"\n✓ Z pierścieniami:    {molecules_with_rings:,} ({molecules_with_rings/len(df)*100:.1f}%)")
print(f"✓ Bez pierścieni:     {molecules_without_rings:,} ({molecules_without_rings/len(df)*100:.1f}%)")

# Wizualizacja
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. Histogram liczby pierścieni
axes[0, 0].bar(ring_dist.index, ring_dist.values, color='steelblue', 
               edgecolor='black', alpha=0.8, linewidth=1)
axes[0, 0].set_xlabel('Liczba pierścieni', fontsize=12, fontweight='bold')
axes[0, 0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0, 0].set_title('Rozkład liczby pierścieni', fontsize=13, fontweight='bold')
axes[0, 0].grid(axis='y', alpha=0.3)

# 2. Pie chart - z/bez pierścieni
axes[0, 1].pie([molecules_with_rings, molecules_without_rings],
               labels=['Z pierścieniami', 'Bez pierścieni'],
               autopct='%1.1f%%', colors=['#3498db', '#95a5a6'],
               startangle=90, textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[0, 1].set_title('Obecność pierścieni', fontsize=13, fontweight='bold')

# 3. Aromatyczne vs Alifatyczne
ring_types = ['aromatic_rings', 'aliphatic_rings', 'saturated_rings']
ring_labels = ['Aromatyczne', 'Alifatyczne', 'Nasycone']
ring_colors = ['#e74c3c', '#2ecc71', '#f39c12']

for ring_type, label, color in zip(ring_types, ring_labels, ring_colors):
    data = df[ring_type].value_counts().sort_index()
    axes[1, 0].plot(data.index, data.values, marker='o', linewidth=2, 
                    label=label, color=color, markersize=8)

axes[1, 0].set_xlabel('Liczba pierścieni', fontsize=12, fontweight='bold')
axes[1, 0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[1, 0].set_title('Typy pierścieni', fontsize=13, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3)

# 4. Heterocykle
hetero_dist = df['hetero_rings'].value_counts().sort_index()
axes[1, 1].bar(hetero_dist.index, hetero_dist.values, color='#9b59b6', 
               edgecolor='black', alpha=0.8, linewidth=1)
axes[1, 1].set_xlabel('Liczba heterocykli', fontsize=12, fontweight='bold')
axes[1, 1].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[1, 1].set_title('Rozkład heterocykli', fontsize=13, fontweight='bold')
axes[1, 1].grid(axis='y', alpha=0.3)

plt.suptitle('Analiza struktur pierścieniowych', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_rings_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_rings_analysis.png")
plt.show()

# ============================================================
# 3. ROTOWALNE WIĄZANIA (FLEXIBILITY)
# ============================================================

print("\n" + "="*70)
print("🔄 ANALIZA ROTOWALNYCH WIĄZAŃ (FLEXIBILITY)")
print("="*70)

def calculate_rotatable_bonds(smiles):
    """Oblicza liczbę rotowalnych wiązań"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    return rdMolDescriptors.CalcNumRotatableBonds(mol)

df['rotatable_bonds'] = df['SMILES'].apply(calculate_rotatable_bonds)

print("\n📊 Statystyki rotowalnych wiązań:")
print(f"  • Średnia:   {df['rotatable_bonds'].mean():.2f} wiązań")
print(f"  • Mediana:   {df['rotatable_bonds'].median():.0f} wiązań")
print(f"  • Max:       {df['rotatable_bonds'].max():.0f} wiązań")
print(f"  • Std:       {df['rotatable_bonds'].std():.2f}")

# Rozkład
print("\n📈 Rozkład rotowalnych wiązań:")
rotatable_dist = df['rotatable_bonds'].value_counts().sort_index()
for count, freq in rotatable_dist.items():
    pct = (freq / len(df)) * 100
    if pct > 1:  # Pokaż tylko istotne
        bar = '█' * int(pct / 2)
        print(f"  {count:2.0f} wiązań: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")

# Klasyfikacja flexibility
rigid = (df['rotatable_bonds'] == 0).sum()
semi_flexible = ((df['rotatable_bonds'] > 0) & (df['rotatable_bonds'] <= 3)).sum()
flexible = (df['rotatable_bonds'] > 3).sum()

print(f"\n📋 Klasyfikacja elastyczności:")
print(f"  • Sztywne (0 wiązań):        {rigid:,} ({rigid/len(df)*100:.1f}%)")
print(f"  • Semi-elastyczne (1-3):     {semi_flexible:,} ({semi_flexible/len(df)*100:.1f}%)")
print(f"  • Elastyczne (4+):           {flexible:,} ({flexible/len(df)*100:.1f}%)")

# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. Histogram
axes[0].hist(df['rotatable_bonds'], bins=range(0, int(df['rotatable_bonds'].max())+2), 
             color='teal', edgecolor='black', alpha=0.7)
axes[0].axvline(df['rotatable_bonds'].mean(), color='red', linestyle='--', 
                linewidth=2, label=f"Średnia: {df['rotatable_bonds'].mean():.2f}")
axes[0].set_xlabel('Liczba rotowalnych wiązań', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład rotowalnych wiązań', fontsize=13, fontweight='bold')
axes[0].legend()
axes[0].grid(axis='y', alpha=0.3)

# 2. Pie chart - klasyfikacja
flexibility_labels = ['Sztywne\n(0)', 'Semi-elastyczne\n(1-3)', 'Elastyczne\n(4+)']
flexibility_values = [rigid, semi_flexible, flexible]
flexibility_colors = ['#e74c3c', '#f39c12', '#2ecc71']

axes[1].pie(flexibility_values, labels=flexibility_labels, autopct='%1.1f%%',
            colors=flexibility_colors, startangle=90,
            textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[1].set_title('Klasyfikacja elastyczności molekuł', fontsize=13, fontweight='bold')

# 3. Zależność: rotatable bonds vs masa molowa
axes[2].scatter(df['MolWt'], df['rotatable_bonds'], alpha=0.3, s=10, color='darkblue')
axes[2].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[2].set_ylabel('Liczba rotowalnych wiązań', fontsize=12, fontweight='bold')
axes[2].set_title('Elastyczność vs Masa molowa', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('qm9_rotatable_bonds.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_rotatable_bonds.png")
plt.show()

# ============================================================
# ============================================================
# 4. FRAKCJA SP3 (SATURACJA VS AROMATYCZNOŚĆ)
# ============================================================

print("\n" + "="*70)
print("🧬 ANALIZA FRAKCJI SP3 (SATURACJA)")
print("="*70)

def calculate_sp3_fraction(smiles):
    """Oblicza frakcję węgli sp3"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    # Liczba węgli sp3
    num_sp3 = 0
    num_carbon = 0

    for atom in mol.GetAtoms():
        if atom.GetSymbol() == 'C':
            num_carbon += 1
            if atom.GetHybridization() == Chem.HybridizationType.SP3:
                num_sp3 += 1

    if num_carbon == 0:
        return 0.0

    return num_sp3 / num_carbon

df['frac_csp3'] = df['SMILES'].apply(calculate_sp3_fraction)

# ============================================================
# 5. CHIRALNOŚĆ
# ============================================================

print("\n" + "="*70)
print("🔀 ANALIZA CHIRALNOŚCI")
print("="*70)

def analyze_chirality(smiles):
    """Analiza centrów chiralnych"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    # Znajdź centra chiralne
    Chem.AssignStereochemistry(mol, cleanIt=True, force=True)
    chiral_centers = Chem.FindMolChiralCenters(mol, includeUnassigned=True)

    return {
        'num_chiral_centers': len(chiral_centers),
        'has_chirality': len(chiral_centers) > 0
    }

print("Analizuję centra chiralne...")
chirality_data = df['SMILES'].apply(analyze_chirality)
df_chirality = chirality_data.apply(pd.Series)

df['num_chiral_centers'] = df_chirality['num_chiral_centers']
df['has_chirality'] = df_chirality['has_chirality']

# Statystyki
chiral_molecules = df['has_chirality'].sum()
achiral_molecules = (~df['has_chirality']).sum()

print(f"\n📊 Statystyki chiralności:")
print(f"  • Cząsteczki chiralne:    {chiral_molecules:,} ({chiral_molecules/len(df)*100:.1f}%)")
print(f"  • Cząsteczki achiralne:   {achiral_molecules:,} ({achiral_molecules/len(df)*100:.1f}%)")

print(f"\n📈 Rozkład liczby centrów chiralnych:")
chiral_dist = df['num_chiral_centers'].value_counts().sort_index()
for count, freq in chiral_dist.items():
    pct = (freq / len(df)) * 100
    if pct > 0.5:
        bar = '█' * int(pct / 2)
        print(f"  {count:1.0f} centrów: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")

# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. Pie chart - chiralne vs achiralne
axes[0].pie([chiral_molecules, achiral_molecules],
            labels=['Chiralne', 'Achiralne'],
            autopct='%1.1f%%', colors=['#e74c3c', '#95a5a6'],
            startangle=90, textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[0].set_title('Obecność chiralności', fontsize=13, fontweight='bold')

# 2. Histogram liczby centrów
chiral_only = df[df['has_chirality']]
if len(chiral_only) > 0:
    axes[1].hist(chiral_only['num_chiral_centers'], 
                 bins=range(0, int(chiral_only['num_chiral_centers'].max())+2),
                 color='coral', edgecolor='black', alpha=0.7)
    axes[1].set_xlabel('Liczba centrów chiralnych', fontsize=12, fontweight='bold')
    axes[1].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
    axes[1].set_title('Rozkład centrów chiralnych\n(tylko molekuły chiralne)', 
                      fontsize=13, fontweight='bold')
    axes[1].grid(axis='y', alpha=0.3)

# 3. Chiralność vs złożoność
axes[2].scatter(df['HeavyAtoms'], df['num_chiral_centers'], alpha=0.3, s=10, color='darkred')
axes[2].set_xlabel('Liczba atomów ciężkich', fontsize=12, fontweight='bold')
axes[2].set_ylabel('Liczba centrów chiralnych', fontsize=12, fontweight='bold')
axes[2].set_title('Chiralność vs Złożoność molekuły', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('qm9_chirality.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_chirality.png")
plt.show()

# ============================================================
# 6. PODSUMOWANIE STRUKTURALNE
# ============================================================

print("\n" + "="*70)
print("📋 PODSUMOWANIE ANALIZY STRUKTURALNEJ")
print("="*70)

print(f"""
🔗 WIĄZANIA:
  • Single:     Średnio {df['single_bonds'].mean():.2f} wiązań
  • Double:     Średnio {df['double_bonds'].mean():.2f} wiązań
  • Triple:     Średnio {df['triple_bonds'].mean():.2f} wiązań
  • Aromatic:   Średnio {df['aromatic_bonds'].mean():.2f} wiązań

⭕ PIERŚCIENIE:
  • Średnia liczba pierścieni:        {df['total_rings'].mean():.2f}
  • Molekuły z pierścieniami:          {molecules_with_rings:,} ({molecules_with_rings/len(df)*100:.1f}%)
  • Średnia pierścieni aromatycznych:  {df['aromatic_rings'].mean():.2f}
  • Średnia heterocykli:               {df['hetero_rings'].mean():.2f}

🔄 FLEXIBILITY:
  • Średnia rotowalnych wiązań:  {df['rotatable_bonds'].mean():.2f}
  • Molekuły sztywne (0):        {rigid:,} ({rigid/len(df)*100:.1f}%)
  • Molekuły elastyczne (4+):    {flexible:,} ({flexible/len(df)*100:.1f}%)

🧬 NASYCENIE:
  • Średnia frakcja sp³:         {df['frac_csp3'].mean():.3f}
  • Min frakcja sp³:             {df['frac_csp3'].min():.3f}
  • Max frakcja sp³:             {df['frac_csp3'].max():.3f}

🔀 CHIRALNOŚĆ:
  • Molekuły chiralne:           {chiral_molecules:,} ({chiral_molecules/len(df)*100:.1f}%)
  • Średnia centrów chiralnych:  {df['num_chiral_centers'].mean():.2f}

✅ CHARAKTERYSTYKA:
  • QM9 zawiera głównie małe, względnie sztywne molekuły
  • Dominują struktury z 0-2 pierścieniami
  • Mix molekuł aromatycznych i alifatycznych
  • Niska chiralność (typowe dla małych molekuł)
""")

print("="*70)
print("✅ ANALIZA STRUKTURALNA ZAKOŃCZONA!")
print("="*70)

======================================================================
ANALIZA STRUKTURALNA - QM9 DATASET
======================================================================

🔗 ETAP 1: ANALIZA TYPÓW WIĄZAŃ
----------------------------------------------------------------------
Analizuję typy wiązań w cząsteczkach...

📊 Statystyki typów wiązań:
            single       double       triple     aromatic
count  1000.000000  1000.000000  1000.000000  1000.000000
mean      7.502000     0.613000     0.268000     0.974000
std       2.624287     0.749528     0.548157     2.129756
min       0.000000     0.000000     0.000000     0.000000
25%       6.000000     0.000000     0.000000     0.000000
50%       8.000000     0.000000     0.000000     0.000000
75%       9.000000     1.000000     0.000000     0.000000
max      13.000000     3.000000     4.000000    10.000000

📈 Średnia liczba wiązań każdego typu:
  • single    : Średnio 7.50 | Molekuły z tym wiązaniem: 997 (99.7%)
  • double    : Średnio 0.61 | Molekuły z tym wiązaniem: 472 (47.2%)
  • triple    : Średnio 0.27 | Molekuły z tym wiązaniem: 222 (22.2%)
  • aromatic  : Średnio 0.97 | Molekuły z tym wiązaniem: 179 (17.9%)

✅ Zapisano: qm9_bond_types.png

======================================================================
⭕ ANALIZA PIERŚCIENI
======================================================================
Analizuję struktury pierścieniowe...

📊 Statystyki pierścieni:
       total_rings  aromatic_rings  aliphatic_rings  saturated_rings  \
count  1000.000000     1000.000000      1000.000000      1000.000000   
mean      1.696000        0.189000         1.507000         1.333000   
std       1.153653        0.416476         1.282028         1.302231   
min       0.000000        0.000000         0.000000         0.000000   
25%       1.000000        0.000000         1.000000         0.000000   
50%       1.000000        0.000000         1.000000         1.000000   
75%       2.000000        0.000000         2.000000         2.000000   
max       7.000000        2.000000         7.000000         7.000000   

       hetero_rings  aromatic_hetero_rings  
count   1000.000000            1000.000000  
mean       1.018000               0.188000  
std        0.893574               0.415727  
min        0.000000               0.000000  
25%        0.000000               0.000000  
50%        1.000000               0.000000  
75%        1.000000               0.000000  
max        6.000000               2.000000  

📈 Rozkład liczby pierścieni:
  0 pierścieni:     98 cząsteczek ( 9.80%) ███
  1 pierścieni:    420 cząsteczek (42.00%) ██████████████
  2 pierścieni:    271 cząsteczek (27.10%) █████████
  3 pierścieni:    146 cząsteczek (14.60%) ████
  4 pierścieni:     38 cząsteczek ( 3.80%) █
  5 pierścieni:     21 cząsteczek ( 2.10%) 
  6 pierścieni:      3 cząsteczek ( 0.30%) 
  7 pierścieni:      3 cząsteczek ( 0.30%) 

✓ Z pierścieniami:    902 (90.2%)
✓ Bez pierścieni:     98 (9.8%)

✅ Zapisano: qm9_rings_analysis.png

======================================================================
🔄 ANALIZA ROTOWALNYCH WIĄZAŃ (FLEXIBILITY)
======================================================================

📊 Statystyki rotowalnych wiązań:
  • Średnia:   0.89 wiązań
  • Mediana:   1 wiązań
  • Max:       5 wiązań
  • Std:       1.07

📈 Rozkład rotowalnych wiązań:
   0 wiązań:    475 cząsteczek (47.50%) ███████████████████████
   1 wiązań:    281 cząsteczek (28.10%) ██████████████
   2 wiązań:    159 cząsteczek (15.90%) ███████
   3 wiązań:     55 cząsteczek ( 5.50%) ██
   4 wiązań:     24 cząsteczek ( 2.40%) █

📋 Klasyfikacja elastyczności:
  • Sztywne (0 wiązań):        475 (47.5%)
  • Semi-elastyczne (1-3):     495 (49.5%)
  • Elastyczne (4+):           30 (3.0%)

✅ Zapisano: qm9_rotatable_bonds.png

======================================================================
🧬 ANALIZA FRAKCJI SP3 (SATURACJA)
======================================================================

======================================================================
🔀 ANALIZA CHIRALNOŚCI
======================================================================
Analizuję centra chiralne...

📊 Statystyki chiralności:
  • Cząsteczki chiralne:    696 (69.6%)
  • Cząsteczki achiralne:   304 (30.4%)

📈 Rozkład liczby centrów chiralnych:
  0 centrów:    304 cząsteczek (30.40%) ███████████████
  1 centrów:    184 cząsteczek (18.40%) █████████
  2 centrów:    203 cząsteczek (20.30%) ██████████
  3 centrów:    133 cząsteczek (13.30%) ██████
  4 centrów:    111 cząsteczek (11.10%) █████
  5 centrów:     44 cząsteczek ( 4.40%) ██
  6 centrów:     17 cząsteczek ( 1.70%) 

✅ Zapisano: qm9_chirality.png

======================================================================
📋 PODSUMOWANIE ANALIZY STRUKTURALNEJ
======================================================================

🔗 WIĄZANIA:
  • Single:     Średnio 7.50 wiązań
  • Double:     Średnio 0.61 wiązań
  • Triple:     Średnio 0.27 wiązań
  • Aromatic:   Średnio 0.97 wiązań

⭕ PIERŚCIENIE:
  • Średnia liczba pierścieni:        1.70
  • Molekuły z pierścieniami:          902 (90.2%)
  • Średnia pierścieni aromatycznych:  0.19
  • Średnia heterocykli:               1.02

🔄 FLEXIBILITY:
  • Średnia rotowalnych wiązań:  0.89
  • Molekuły sztywne (0):        475 (47.5%)
  • Molekuły elastyczne (4+):    30 (3.0%)

🧬 NASYCENIE:
  • Średnia frakcja sp³:         0.695
  • Min frakcja sp³:             0.000
  • Max frakcja sp³:             1.000

🔀 CHIRALNOŚĆ:
  • Molekuły chiralne:           696 (69.6%)
  • Średnia centrów chiralnych:  1.78

✅ CHARAKTERYSTYKA:
  • QM9 zawiera głównie małe, względnie sztywne molekuły
  • Dominują struktury z 0-2 pierścieniami
  • Mix molekuł aromatycznych i alifatycznych
  • Niska chiralność (typowe dla małych molekuł)

======================================================================
✅ ANALIZA STRUKTURALNA ZAKOŃCZONA!
======================================================================

5. WŁAŚCIWOŚCI FIZYKOCHEMICZNE (RDKit) ⚗️¶

• LogP - lipofilowość, rozkład, outliers
• TPSA - polar surface area, drug-likeness
• HBD/HBA - donory i akceptory wodoru
• Molar Refractivity - objętość molekularna
• Reguła Lipińskiego - ile struktur spełnia Ro5
• Korelacje między deskryptorami - heatmap

In [26]:

# ============================================================
# KOMÓRKA 11: WŁAŚCIWOŚCI FIZYKOCHEMICZNE (RDKit)
# ============================================================

print("="*70)
print("WŁAŚCIWOŚCI FIZYKOCHEMICZNE - QM9 DATASET")
print("="*70)

# ============================================================
# 1. OBLICZANIE DESKRYPTORÓW RDKIT
# ============================================================

print("\n⚗️ ETAP 1: OBLICZANIE DESKRYPTORÓW FIZYKOCHEMICZNYCH")
print("-" * 70)

from rdkit import Chem
from rdkit.Chem import Descriptors, Crippen, Lipinski

def calculate_physicochemical_properties(smiles):
    """Oblicza wszystkie kluczowe właściwości fizykochemiczne"""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    try:
        return {
            'LogP': Crippen.MolLogP(mol),
            'MolMR': Crippen.MolMR(mol),  # Molar Refractivity
            'TPSA': Descriptors.TPSA(mol),
            'HBD': Lipinski.NumHDonors(mol),  # H-Bond Donors
            'HBA': Lipinski.NumHAcceptors(mol),  # H-Bond Acceptors
            'LabuteASA': Descriptors.LabuteASA(mol),  # Surface Area
            'BalabanJ': Descriptors.BalabanJ(mol),  # Topological index
        }
    except:
        return None

print("Obliczam właściwości fizykochemiczne...")
props_data = df['SMILES'].apply(calculate_physicochemical_properties)
df_props = props_data.apply(pd.Series)

# Dodanie do głównego DataFrame
for col in df_props.columns:
    df[col] = df_props[col]

print(f"\n✅ Obliczono {len(df_props.columns)} właściwości dla {len(df)} molekuł")

print("\n📊 Statystyki właściwości fizykochemicznych:")
print(df_props.describe())

# ============================================================
# 2. ANALIZA LogP (LIPOFILOWOŚĆ)
# ============================================================

print("\n" + "="*70)
print("💧 ANALIZA LogP (LIPOFILOWOŚĆ)")
print("="*70)

print("\n📊 Statystyki LogP:")
print(f"  • Średnia:    {df['LogP'].mean():.3f}")
print(f"  • Mediana:    {df['LogP'].median():.3f}")
print(f"  • Std:        {df['LogP'].std():.3f}")
print(f"  • Min:        {df['LogP'].min():.3f}")
print(f"  • Max:        {df['LogP'].max():.3f}")

# Percentyle
print("\n📈 Percentyle LogP:")
for p in [10, 25, 50, 75, 90]:
    val = df['LogP'].quantile(p/100)
    print(f"  • {p:2d}%:  {val:.3f}")

# Klasyfikacja lipofilowości
very_hydrophilic = (df['LogP'] < -1).sum()
hydrophilic = ((df['LogP'] >= -1) & (df['LogP'] < 1)).sum()
moderate = ((df['LogP'] >= 1) & (df['LogP'] < 3)).sum()
lipophilic = ((df['LogP'] >= 3) & (df['LogP'] < 5)).sum()
very_lipophilic = (df['LogP'] >= 5).sum()

print(f"\n📋 Klasyfikacja lipofilowości:")
print(f"  • Bardzo hydrofilowe (LogP < -1):    {very_hydrophilic:,} ({very_hydrophilic/len(df)*100:.1f}%)")
print(f"  • Hydrofilowe (-1 ≤ LogP < 1):       {hydrophilic:,} ({hydrophilic/len(df)*100:.1f}%)")
print(f"  • Umiarkowane (1 ≤ LogP < 3):        {moderate:,} ({moderate/len(df)*100:.1f}%)")
print(f"  • Lipofilowe (3 ≤ LogP < 5):         {lipophilic:,} ({lipophilic/len(df)*100:.1f}%)")
print(f"  • Bardzo lipofilowe (LogP ≥ 5):      {very_lipophilic:,} ({very_lipophilic/len(df)*100:.1f}%)")

# Outliers (IQR method)
Q1 = df['LogP'].quantile(0.25)
Q3 = df['LogP'].quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
outliers_logp = df[(df['LogP'] < lower_bound) | (df['LogP'] > upper_bound)]

print(f"\n⚠️  Outliers LogP (metoda IQR):")
print(f"  • Liczba outliers: {len(outliers_logp)} ({len(outliers_logp)/len(df)*100:.2f}%)")
print(f"  • Zakres normalny: [{lower_bound:.2f}, {upper_bound:.2f}]")

# Wizualizacja
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. Histogram
axes[0, 0].hist(df['LogP'], bins=60, color='dodgerblue', edgecolor='black', alpha=0.7)
axes[0, 0].axvline(df['LogP'].mean(), color='red', linestyle='--', linewidth=2, 
                   label=f"Średnia: {df['LogP'].mean():.2f}")
axes[0, 0].axvline(df['LogP'].median(), color='green', linestyle='--', linewidth=2,
                   label=f"Mediana: {df['LogP'].median():.2f}")
axes[0, 0].set_xlabel('LogP', fontsize=12, fontweight='bold')
axes[0, 0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[0, 0].set_title('Rozkład LogP', fontsize=13, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(axis='y', alpha=0.3)

# 2. Boxplot z outliers
bp = axes[0, 1].boxplot([df['LogP'].dropna()], labels=['LogP'], 
                        patch_artist=True, widths=0.5, showfliers=True)
bp['boxes'][0].set_facecolor('lightblue')
axes[0, 1].set_ylabel('LogP', fontsize=12, fontweight='bold')
axes[0, 1].set_title('Boxplot LogP (z outliers)', fontsize=13, fontweight='bold')
axes[0, 1].grid(axis='y', alpha=0.3)

# 3. Klasyfikacja - barplot
categories = ['Bardzo\nhydrofilowe', 'Hydrofilowe', 'Umiarkowane', 
              'Lipofilowe', 'Bardzo\nlipofilowe']
values = [very_hydrophilic, hydrophilic, moderate, lipophilic, very_lipophilic]
colors_logp = ['#3498db', '#5dade2', '#aed6f1', '#f8b739', '#e67e22']

bars = axes[1, 0].bar(range(len(categories)), values, color=colors_logp, 
                      edgecolor='black', alpha=0.8, linewidth=1)
axes[1, 0].set_xticks(range(len(categories)))
axes[1, 0].set_xticklabels(categories, fontsize=9)
axes[1, 0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[1, 0].set_title('Klasyfikacja lipofilowości', fontsize=13, fontweight='bold')
axes[1, 0].grid(axis='y', alpha=0.3)

for bar, val in zip(bars, values):
    axes[1, 0].text(bar.get_x() + bar.get_width()/2, val + 10, 
                    f'{val}', ha='center', fontsize=9, fontweight='bold')

# 4. LogP vs Masa molowa
axes[1, 1].scatter(df['MolWt'], df['LogP'], alpha=0.3, s=10, color='darkblue')
axes[1, 1].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[1, 1].set_ylabel('LogP', fontsize=12, fontweight='bold')
axes[1, 1].set_title('LogP vs Masa molowa', fontsize=13, fontweight='bold')
axes[1, 1].grid(alpha=0.3)

plt.suptitle('Analiza lipofilowości (LogP)', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_logp_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_logp_analysis.png")
plt.show()

# ============================================================
# 3. ANALIZA TPSA (POLAR SURFACE AREA)
# ============================================================

print("\n" + "="*70)
print("🧊 ANALIZA TPSA (TOPOLOGICAL POLAR SURFACE AREA)")
print("="*70)

print("\n📊 Statystyki TPSA:")
print(f"  • Średnia:    {df['TPSA'].mean():.2f} Ų")
print(f"  • Mediana:    {df['TPSA'].median():.2f} Ų")
print(f"  • Std:        {df['TPSA'].std():.2f} Ų")
print(f"  • Max:        {df['TPSA'].max():.2f} Ų")

# Drug-likeness według TPSA
# TPSA < 140 Ų - dobra przepuszczalność BBB (blood-brain barrier)
# TPSA < 90 Ų - optymalna dla leków CNS
bbb_permeable = (df['TPSA'] < 140).sum()
cns_optimal = (df['TPSA'] < 90).sum()

print(f"\n💊 Drug-likeness (TPSA):")
print(f"  • BBB permeable (< 140 Ų):    {bbb_permeable:,} ({bbb_permeable/len(df)*100:.1f}%)")
print(f"  • CNS optimal (< 90 Ų):       {cns_optimal:,} ({cns_optimal/len(df)*100:.1f}%)")

# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. Histogram
axes[0].hist(df['TPSA'], bins=50, color='coral', edgecolor='black', alpha=0.7)
axes[0].axvline(90, color='green', linestyle='--', linewidth=2, label='CNS optimal (90)')
axes[0].axvline(140, color='orange', linestyle='--', linewidth=2, label='BBB limit (140)')
axes[0].set_xlabel('TPSA (Ų)', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład TPSA', fontsize=13, fontweight='bold')
axes[0].legend()
axes[0].grid(axis='y', alpha=0.3)

# 2. TPSA vs LogP
axes[1].scatter(df['LogP'], df['TPSA'], alpha=0.3, s=10, color='purple')
axes[1].set_xlabel('LogP', fontsize=12, fontweight='bold')
axes[1].set_ylabel('TPSA (Ų)', fontsize=12, fontweight='bold')
axes[1].set_title('TPSA vs LogP', fontsize=13, fontweight='bold')
axes[1].grid(alpha=0.3)

# 3. Pie chart - drug-likeness
categories_tpsa = ['CNS optimal\n(<90)', 'BBB permeable\n(90-140)', 'High TPSA\n(>140)']
values_tpsa = [cns_optimal, bbb_permeable - cns_optimal, len(df) - bbb_permeable]
colors_tpsa = ['#2ecc71', '#f39c12', '#e74c3c']

axes[2].pie(values_tpsa, labels=categories_tpsa, autopct='%1.1f%%',
            colors=colors_tpsa, startangle=90, 
            textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[2].set_title('TPSA Drug-likeness', fontsize=13, fontweight='bold')

plt.tight_layout()
plt.savefig('qm9_tpsa_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_tpsa_analysis.png")
plt.show()

# ============================================================
# 4. WIĄZANIA WODOROWE (HBD/HBA)
# ============================================================

print("\n" + "="*70)
print("🔗 ANALIZA WIĄZAŃ WODOROWYCH (HBD/HBA)")
print("="*70)

print("\n📊 Statystyki:")
print(f"  • Średnia HBD (donory):     {df['HBD'].mean():.2f}")
print(f"  • Średnia HBA (akceptory):  {df['HBA'].mean():.2f}")
print(f"  • Max HBD:                  {df['HBD'].max():.0f}")
print(f"  • Max HBA:                  {df['HBA'].max():.0f}")

# Rozkład
print("\n📈 Rozkład HBD:")
hbd_dist = df['HBD'].value_counts().sort_index()
for count, freq in hbd_dist.items():
    pct = (freq / len(df)) * 100
    if pct > 1:
        bar = '█' * int(pct / 3)
        print(f"  {count:.0f} donorów: {freq:5,} ({pct:5.2f}%) {bar}")

print("\n📈 Rozkład HBA:")
hba_dist = df['HBA'].value_counts().sort_index()
for count, freq in hba_dist.items():
    pct = (freq / len(df)) * 100
    if pct > 1:
        bar = '█' * int(pct / 3)
        print(f"  {count:.0f} akceptorów: {freq:5,} ({pct:5.2f}%) {bar}")

# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. HBD histogram
axes[0].bar(hbd_dist.index, hbd_dist.values, color='#3498db', 
            edgecolor='black', alpha=0.8, linewidth=1)
axes[0].set_xlabel('Liczba donorów H (HBD)', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład donorów wodoru', fontsize=13, fontweight='bold')
axes[0].grid(axis='y', alpha=0.3)

# 2. HBA histogram
axes[1].bar(hba_dist.index, hba_dist.values, color='#e74c3c', 
            edgecolor='black', alpha=0.8, linewidth=1)
axes[1].set_xlabel('Liczba akceptorów H (HBA)', fontsize=12, fontweight='bold')
axes[1].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[1].set_title('Rozkład akceptorów wodoru', fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)

# 3. HBD vs HBA scatter
axes[2].scatter(df['HBD'], df['HBA'], alpha=0.4, s=15, color='green', edgecolors='black', linewidth=0.3)
axes[2].set_xlabel('HBD (donory)', fontsize=12, fontweight='bold')
axes[2].set_ylabel('HBA (akceptory)', fontsize=12, fontweight='bold')
axes[2].set_title('HBD vs HBA', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('qm9_hbond_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_hbond_analysis.png")
plt.show()

# ============================================================
# 5. REGUŁA LIPIŃSKIEGO (RULE OF 5)
# ============================================================

print("\n" + "="*70)
print("💊 REGUŁA LIPIŃSKIEGO (RULE OF 5)")
print("="*70)

# Kryteria Lipińskiego:
# 1. Masa molowa ≤ 500 Da
# 2. LogP ≤ 5
# 3. HBD ≤ 5
# 4. HBA ≤ 10

df['Lipinski_MW'] = df['MolWt'] <= 500
df['Lipinski_LogP'] = df['LogP'] <= 5
df['Lipinski_HBD'] = df['HBD'] <= 5
df['Lipinski_HBA'] = df['HBA'] <= 10

df['Lipinski_violations'] = (~df['Lipinski_MW']).astype(int) + \
                             (~df['Lipinski_LogP']).astype(int) + \
                             (~df['Lipinski_HBD']).astype(int) + \
                             (~df['Lipinski_HBA']).astype(int)

df['Lipinski_compliant'] = df['Lipinski_violations'] <= 1  # 0 lub 1 naruszenie = OK

# Statystyki
compliant = df['Lipinski_compliant'].sum()
non_compliant = (~df['Lipinski_compliant']).sum()

print(f"\n📊 Zgodność z regułą Lipińskiego:")
print(f"  • Zgodne (≤1 naruszenie):    {compliant:,} ({compliant/len(df)*100:.1f}%)")
print(f"  • Niezgodne (>1 naruszenie): {non_compliant:,} ({non_compliant/len(df)*100:.1f}%)")

print(f"\n📈 Rozkład naruszeń:")
violations_dist = df['Lipinski_violations'].value_counts().sort_index()
for count, freq in violations_dist.items():
    pct = (freq / len(df)) * 100
    bar = '█' * int(pct / 2)
    print(f"  {count:.0f} naruszeń: {freq:5,} ({pct:5.2f}%) {bar}")

print(f"\n📋 Szczegółowe naruszenia:")
print(f"  • MW > 500:    {(~df['Lipinski_MW']).sum():,}")
print(f"  • LogP > 5:    {(~df['Lipinski_LogP']).sum():,}")
print(f"  • HBD > 5:     {(~df['Lipinski_HBD']).sum():,}")
print(f"  • HBA > 10:    {(~df['Lipinski_HBA']).sum():,}")

# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. Pie chart - compliance
axes[0].pie([compliant, non_compliant],
            labels=['Zgodne\n(Ro5)', 'Niezgodne'],
            autopct='%1.1f%%', colors=['#2ecc71', '#e74c3c'],
            startangle=90, textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[0].set_title('Zgodność z regułą Lipińskiego', fontsize=13, fontweight='bold')

# 2. Barplot naruszeń
axes[1].bar(violations_dist.index, violations_dist.values, color='#9b59b6', 
            edgecolor='black', alpha=0.8, linewidth=1)
axes[1].set_xlabel('Liczba naruszeń Ro5', fontsize=12, fontweight='bold')
axes[1].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[1].set_title('Rozkład naruszeń Ro5', fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)

# 3. Szczegółowe naruszenia
violation_types = ['MW > 500', 'LogP > 5', 'HBD > 5', 'HBA > 10']
violation_counts = [
    (~df['Lipinski_MW']).sum(),
    (~df['Lipinski_LogP']).sum(),
    (~df['Lipinski_HBD']).sum(),
    (~df['Lipinski_HBA']).sum()
]

axes[2].barh(range(len(violation_types)), violation_counts, color='#e74c3c', 
             edgecolor='black', alpha=0.8, linewidth=1)
axes[2].set_yticks(range(len(violation_types)))
axes[2].set_yticklabels(violation_types)
axes[2].set_xlabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[2].set_title('Typy naruszeń Ro5', fontsize=13, fontweight='bold')
axes[2].grid(axis='x', alpha=0.3)

for i, val in enumerate(violation_counts):
    axes[2].text(val + 5, i, f'{val}', va='center', fontweight='bold')

plt.tight_layout()
plt.savefig('qm9_lipinski_ro5.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_lipinski_ro5.png")
plt.show()

# ============================================================
# 6. KORELACJE MIĘDZY DESKRYPTORAMI
# ============================================================

print("\n" + "="*70)
print("🔗 KORELACJE MIĘDZY DESKRYPTORAMI")
print("="*70)

# Wybór deskryptorów do analizy
descriptors_to_correlate = ['MolWt', 'LogP', 'TPSA', 'HBD', 'HBA', 'MolMR', 
                            'HeavyAtoms', 'rotatable_bonds', 'frac_csp3']

df_corr = df[descriptors_to_correlate].copy()

# Macierz korelacji
corr_matrix = df_corr.corr()

print("\n📊 Najsilniejsze korelacje:")
# Ekstakcja korelacji (bez diagonali)
correlations_list = []
for i in range(len(corr_matrix.columns)):
    for j in range(i+1, len(corr_matrix.columns)):
        correlations_list.append({
            'Prop1': corr_matrix.columns[i],
            'Prop2': corr_matrix.columns[j],
            'Correlation': corr_matrix.iloc[i, j]
        })

df_corr_list = pd.DataFrame(correlations_list)
df_corr_list = df_corr_list.sort_values('Correlation', key=abs, ascending=False)

print("\nTop 10 najsilniejszych korelacji:")
for idx, row in df_corr_list.head(10).iterrows():
    print(f"  {row['Prop1']:20s} ↔ {row['Prop2']:20s}  r = {row['Correlation']:+.3f}")

# Wizualizacja
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# 1. Heatmap
import seaborn as sns
sns.heatmap(corr_matrix, annot=True, fmt='.2f', cmap='RdBu_r', 
            center=0, square=True, linewidths=0.5, cbar_kws={"shrink": 0.8},
            vmin=-1, vmax=1, ax=axes[0])
axes[0].set_title('Macierz korelacji deskryptorów', fontsize=13, fontweight='bold')

# 2. Top korelacje - barplot
top10 = df_corr_list.head(10).copy()
top10['Pair'] = top10['Prop1'] + ' ↔\n' + top10['Prop2']
colors_corr = ['#d73027' if x < 0 else '#4575b4' for x in top10['Correlation']]

axes[1].barh(range(len(top10)), top10['Correlation'], color=colors_corr, 
             edgecolor='black', alpha=0.8, linewidth=0.7)
axes[1].set_yticks(range(len(top10)))
axes[1].set_yticklabels(top10['Pair'], fontsize=8)
axes[1].set_xlabel('Współczynnik korelacji', fontsize=12, fontweight='bold')
axes[1].set_title('Top 10 najsilniejszych korelacji', fontsize=13, fontweight='bold')
axes[1].axvline(0, color='black', linewidth=0.8)
axes[1].grid(axis='x', alpha=0.3)
axes[1].invert_yaxis()

for i, val in enumerate(top10['Correlation']):
    axes[1].text(val + 0.02 if val > 0 else val - 0.02, i, 
                 f'{val:.2f}', va='center', ha='left' if val > 0 else 'right',
                 fontweight='bold', fontsize=8)

plt.tight_layout()
plt.savefig('qm9_descriptors_correlation.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_descriptors_correlation.png")
plt.show()

# ============================================================
# 7. PODSUMOWANIE
# ============================================================

print("\n" + "="*70)
print("📋 PODSUMOWANIE WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH")
print("="*70)

print(f"""
💧 LIPOFILOWOŚĆ (LogP):
  • Średnia:              {df['LogP'].mean():.3f}
  • Zakres:               {df['LogP'].min():.2f} - {df['LogP'].max():.2f}
  • Hydrofilowe (<1):     {(df['LogP'] < 1).sum():,} ({(df['LogP'] < 1).sum()/len(df)*100:.1f}%)
  • Lipofilowe (≥3):      {(df['LogP'] >= 3).sum():,} ({(df['LogP'] >= 3).sum()/len(df)*100:.1f}%)

🧊 TPSA (Polar Surface Area):
  • Średnia:              {df['TPSA'].mean():.2f} Ų
  • BBB permeable:        {bbb_permeable:,} ({bbb_permeable/len(df)*100:.1f}%)
  • CNS optimal:          {cns_optimal:,} ({cns_optimal/len(df)*100:.1f}%)

🔗 WIĄZANIA WODOROWE:
  • Średnia HBD:          {df['HBD'].mean():.2f}
  • Średnia HBA:          {df['HBA'].mean():.2f}
  • Max HBD:              {df['HBD'].max():.0f}
  • Max HBA:              {df['HBA'].max():.0f}

💊 REGUŁA LIPIŃSKIEGO (Ro5):
  • Zgodne molekuły:      {compliant:,} ({compliant/len(df)*100:.1f}%)
  • 0 naruszeń:           {(df['Lipinski_violations']==0).sum():,}
  • 1 naruszenie:         {(df['Lipinski_violations']==1).sum():,}
  • >1 naruszeń:          {(df['Lipinski_violations']>1).sum():,}

🔗 KORELACJE:
  • Najsilniejsza:        {df_corr_list.iloc[0]['Prop1']} ↔ {df_corr_list.iloc[0]['Prop2']} (r={df_corr_list.iloc[0]['Correlation']:.3f})
  • Silne korelacje (|r|>0.7): {(df_corr_list['Correlation'].abs() > 0.7).sum()}

✅ WNIOSKI:
  • QM9 zawiera głównie małe, relatywnie hydrofilowe molekuły
  • Wysoka zgodność z regułą Lipińskiego ({compliant/len(df)*100:.0f}%)
  • Niska TPSA → dobra przepuszczalność błon
  • Deskryptory silnie skorelowane z masą molową
""")

print("="*70)
print("✅ ANALIZA WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH ZAKOŃCZONA!")
print("="*70)

======================================================================
WŁAŚCIWOŚCI FIZYKOCHEMICZNE - QM9 DATASET
======================================================================

⚗️ ETAP 1: OBLICZANIE DESKRYPTORÓW FIZYKOCHEMICZNYCH
----------------------------------------------------------------------
Obliczam właściwości fizykochemiczne...

✅ Obliczono 7 właściwości dla 1000 molekuł

📊 Statystyki właściwości fizykochemicznych:
              LogP        MolMR         TPSA          HBD          HBA  \
count  1000.000000  1000.000000  1000.000000  1000.000000  1000.000000   
mean      0.104381    31.821830    37.802210     0.972000     2.053000   
std       1.173836     3.606108    21.399043     0.889948     1.093337   
min      -3.467400    14.868400     0.000000     0.000000     0.000000   
25%      -0.660365    29.675525    21.510000     0.000000     1.000000   
50%       0.158600    31.994000    36.510000     1.000000     2.000000   
75%       0.923925    34.185250    52.512500     2.000000     3.000000   
max       3.222700    41.365000   112.360000     4.000000     5.000000   

         LabuteASA     BalabanJ  
count  1000.000000  1000.000000  
mean     52.579936     2.448441  
std       3.722869     0.525115  
min      24.605606     0.000000  
25%      51.782555     2.071725  
50%      53.587200     2.269522  
75%      54.820772     2.693324  
max      58.088080     4.837063  

======================================================================
💧 ANALIZA LogP (LIPOFILOWOŚĆ)
======================================================================

📊 Statystyki LogP:
  • Średnia:    0.104
  • Mediana:    0.159
  • Std:        1.174
  • Min:        -3.467
  • Max:        3.223

📈 Percentyle LogP:
  • 10%:  -1.387
  • 25%:  -0.660
  • 50%:  0.159
  • 75%:  0.924
  • 90%:  1.561

📋 Klasyfikacja lipofilowości:
  • Bardzo hydrofilowe (LogP < -1):    179 (17.9%)
  • Hydrofilowe (-1 ≤ LogP < 1):       592 (59.2%)
  • Umiarkowane (1 ≤ LogP < 3):        227 (22.7%)
  • Lipofilowe (3 ≤ LogP < 5):         2 (0.2%)
  • Bardzo lipofilowe (LogP ≥ 5):      0 (0.0%)

⚠️  Outliers LogP (metoda IQR):
  • Liczba outliers: 7 (0.70%)
  • Zakres normalny: [-3.04, 3.30]

C:\Users\slast\AppData\Local\Temp\ipykernel_6436\1039049785.py:114: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11.
  bp = axes[0, 1].boxplot([df['LogP'].dropna()], labels=['LogP'],

✅ Zapisano: qm9_logp_analysis.png

======================================================================
🧊 ANALIZA TPSA (TOPOLOGICAL POLAR SURFACE AREA)
======================================================================

📊 Statystyki TPSA:
  • Średnia:    37.80 Ų
  • Mediana:    36.51 Ų
  • Std:        21.40 Ų
  • Max:        112.36 Ų

💊 Drug-likeness (TPSA):
  • BBB permeable (< 140 Ų):    1,000 (100.0%)
  • CNS optimal (< 90 Ų):       991 (99.1%)

✅ Zapisano: qm9_tpsa_analysis.png

======================================================================
🔗 ANALIZA WIĄZAŃ WODOROWYCH (HBD/HBA)
======================================================================

📊 Statystyki:
  • Średnia HBD (donory):     0.97
  • Średnia HBA (akceptory):  2.05
  • Max HBD:                  4
  • Max HBA:                  5

📈 Rozkład HBD:
  0 donorów:   348 (34.80%) ███████████
  1 donorów:   392 (39.20%) █████████████
  2 donorów:   202 (20.20%) ██████
  3 donorów:    56 ( 5.60%) █

📈 Rozkład HBA:
  0 akceptorów:    66 ( 6.60%) ██
  1 akceptorów:   249 (24.90%) ████████
  2 akceptorów:   355 (35.50%) ███████████
  3 akceptorów:   247 (24.70%) ████████
  4 akceptorów:    62 ( 6.20%) ██
  5 akceptorów:    21 ( 2.10%) 

✅ Zapisano: qm9_hbond_analysis.png

======================================================================
💊 REGUŁA LIPIŃSKIEGO (RULE OF 5)
======================================================================

📊 Zgodność z regułą Lipińskiego:
  • Zgodne (≤1 naruszenie):    1,000 (100.0%)
  • Niezgodne (>1 naruszenie): 0 (0.0%)

📈 Rozkład naruszeń:
  0 naruszeń: 1,000 (100.00%) ██████████████████████████████████████████████████

📋 Szczegółowe naruszenia:
  • MW > 500:    0
  • LogP > 5:    0
  • HBD > 5:     0
  • HBA > 10:    0

C:\Users\slast\AppData\Local\Temp\ipykernel_6436\1039049785.py:358: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all Axes decorations.
  plt.tight_layout()

✅ Zapisano: qm9_lipinski_ro5.png

======================================================================
🔗 KORELACJE MIĘDZY DESKRYPTORAMI
======================================================================

📊 Najsilniejsze korelacje:

Top 10 najsilniejszych korelacji:
  MolWt                ↔ HeavyAtoms            r = +0.927
  TPSA                 ↔ HBA                   r = +0.821
  LogP                 ↔ HBD                   r = -0.638
  MolMR                ↔ HeavyAtoms            r = +0.578
  TPSA                 ↔ HBD                   r = +0.572
  LogP                 ↔ MolMR                 r = +0.495
  LogP                 ↔ TPSA                  r = -0.476
  MolWt                ↔ MolMR                 r = +0.470
  HBA                  ↔ MolMR                 r = -0.462
  TPSA                 ↔ frac_csp3             r = -0.447

✅ Zapisano: qm9_descriptors_correlation.png

======================================================================
📋 PODSUMOWANIE WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH
======================================================================

💧 LIPOFILOWOŚĆ (LogP):
  • Średnia:              0.104
  • Zakres:               -3.47 - 3.22
  • Hydrofilowe (<1):     771 (77.1%)
  • Lipofilowe (≥3):      2 (0.2%)

🧊 TPSA (Polar Surface Area):
  • Średnia:              37.80 Ų
  • BBB permeable:        1,000 (100.0%)
  • CNS optimal:          991 (99.1%)

🔗 WIĄZANIA WODOROWE:
  • Średnia HBD:          0.97
  • Średnia HBA:          2.05
  • Max HBD:              4
  • Max HBA:              5

💊 REGUŁA LIPIŃSKIEGO (Ro5):
  • Zgodne molekuły:      1,000 (100.0%)
  • 0 naruszeń:           1,000
  • 1 naruszenie:         0
  • >1 naruszeń:          0

🔗 KORELACJE:
  • Najsilniejsza:        MolWt ↔ HeavyAtoms (r=0.927)
  • Silne korelacje (|r|>0.7): 2

✅ WNIOSKI:
  • QM9 zawiera głównie małe, relatywnie hydrofilowe molekuły
  • Wysoka zgodność z regułą Lipińskiego (100%)
  • Niska TPSA → dobra przepuszczalność błon
  • Deskryptory silnie skorelowane z masą molową

======================================================================
✅ ANALIZA WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH ZAKOŃCZONA!
======================================================================